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Prediction  of  the  behavior  of  neutron  waves  in  heterogeneous  media 
is  performed  using  one-group  diffusion  theory  and  age-diffusion  theory.   The 
one-group  theory  is  used  in  two  different^  finite  geometries  which  allow 
reduction  to  a  one -dimensional  problem  and  to  a  two-dimensional  problem. 
These  cases  result  in  diffusion  kernels,  or  Green's  functions,  for  the  two 
finite  configurations. 

The  theory  is  then  extended  to  include  Fermi-age,  or  continuous 
slowing-down,  theory  for  higher  energy  neutrons.   The  approach  is  similar 
to  the  Feinberg-Galanin  heterogeneous  reactor  theory  except  that  it  is  ap- 
plied to  a  finite  geometry  and  includes  time  dependence.   A  finite  diffu- 
sion kernel  is  obtained  which  is  similar  to  the  results  of  the  simpler 
calculations.   In  addition,  a  finite -medium,  Fermi -age  kernel  results  which 
describes  the  behavior  of  the  slowing-down  neutrons  in  the  finite  geometry. 


Vll 


Both  the  one-group  and  the  age-diffusion  developments  are  used  to 
calculate  numerically  sample  configurations  which  are  suitable  for  experi- 
mental verification.   These  include  rectangular  assemblies  of  graphite  and 
heavy  water  which  have  poison  or  fuel  rods  inserted. 

In  order  to  demonstrate  better  the  improvements  of  this  work  over 
the  Feinberg-Galanin  theory^  an  extension  to  critical  assemblies  is  made. 
This  offers  the  possibility  of  doing  criticality  calculations  for  physically 
small  assemblies  which  are  not  amenable  to  homogenization  techniques. 

Using  the  age-diffusion  theory  results^  it  is  also  shown  that  data 
from  neutron  wave  experiments  may  be  processed  in  such  a  way  as  to  give 
experimental  measurements  of  both  the  diffusion  and  the  slowing-down 
kernel. 
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INTRODUCTION 

In  i960  a  program  was  initiated  at  the  University  of  Florida  in  the 
field  of  neutron  wave  propagation.   The  initial  effort  was  in  the  determina- 
tion of  diffusion  and  thermalization  parameters  in  moderating  media.   It  was 

,  ,1 
shown  by  Perez  and  Uhrig  (1)      that  these  parameters  were  related  to  the 

damping  coefficient  and  phase  shift  per  unit  length  of  the  neutron  wave. 
The  work  of  Hartley  (2)  and  Booth  (5)  confirmed  these  predictions^  and 
lately  Perez  and  Booth  (h)   were  able  to  perform  accixrate  measiurements  of 
thermalization  parameters  in  graphite  using  the  neutron  wave  technique. 
The  above  mentioned  studies  were  performed  in  homogeneous  neutron  moderat- 
ing media.   Some  work  has  been  done  by  Booth  (5.)^  in  heterogeneous  media^ 
and  by  Denning  e_t  al.  (6)  in  two-region  moderating  media.  Denning  et  al. 
were  able  to  show  a  definite  reflected  wave  from  the  interface  between  the 
two  media  which  interferes  with  the  incident  wave. 

The  goal  of  the  present  work  is  to  predict,  as  accurately  as  is  fea- 
sible, the  propagation  of  neutron  waves  in  heterogeneous  multiplying  and 
nonmul tip lying  media.   The  incentive  to  perform  this  work  is  twofold.   From 
the  purely  academic  viewpoint,  the  propagation  of  waves  through  periodic 
structures  is  of  sufficient  interest  in  the  fields  of  quantum  mechanics, 
sound  theory,  solid-state  physics,  and  others  that  extension  of  the 


■""Underlined  numbers  in  parentheses  refer  to  corresponding  numbers 
in  the  list  of  references. 


previously  mentioned  work  in  heterogeneous  media  becomes  desirable.   That 
neutron  wave  propagation  theory  is  applicable  to  these  other  fields  may 
perhaps  seem  reasonable  by  pointing  out  that  the  neutron  wave  propagation 
we  are  concerned  with  here  is  a  macroscopic  phenomenon,  qualitatively 
similar  to  soimd  wave  propagation.   One  simply  replaces  gas  density  with 
neutron  density.   Quantitatively,  the  neutron  wave  propagation  is  the  more 
difficult  problem  since  nuclear  reactors  are  highly  dispersive  and  absorp- 
tive media  in  which  disturbances  of  the  neutron  field  are  qiJickly  damped. 

Secondly,  there  is  a  strong  incentive  from  the  practical  viewpoint 
of  gaining  understanding  about  heterogeneous  reactors.   Traditionally, 
homogenization  techniques  based  on  the  Wigner-Seitz  unit-cell  theory  have 
been  used  to  design  heterogeneous  reactors,  and  it  is  only  recently  that 
theories  accounting  for  the  interactions  between  fuel  plates  are  evolving 
from  the  works  of  Galanin  (7);,  Feinberg  (8),  and  others. 

The  study  of  neutron  wave  propagation  through  heterogeneous  media 
involves  analytical  techniques  and  ass-umptions  similar  to  the  developments 
used  in  heterogeneous  reactor  theory.  Hence  studies  of  this  type  offer  the 
possibilities  of  experimental  tests  of  the  various  theories  in  a  way  which 
enlarges  the  realm  of  application  of  purely  static  experiments.   Also,  any 
disturbance  of  the  neutron  field  in  a  heterogeneous  reactor  can  be  analyzed 
as  a  superposition  of  neutron  waves,  thus  yielding  information  on  possible 
short-range  instabilities  for  a  particular  reactor  design. 

Theoretical  studies  in  this  field  are  q\aite  complicated  because  of 
the  interplay  of  the  effects  of  the  spatial  heterogeneities  produced  by  the 
fuel  plates  and  the  energy  distribution  of  the  neutron  population.   Before 


meaningful  but  costly  experiments  can  be  performed^  the  theoretical  founda- 
tions have  to  be  developed^  which  is  the  goal  set  forth  in  this  disserta- 
tion. 

At  first  sight^  the  obvious  approach  to  the  theoretical  work  is 
simply  to  modify  the  Feinberg-Galanin  theory  to  include  time  dependence  so 
that  it  can  be  applied  to  subcritical  assemblies.   In  practice,  this  is  not 
satisfactory  since  most  reactors  are  rather  large  assemblies.  Because  of 
this  the  present  heterogeneous  reactor  theory  assumes  that  the  moderating 
medium  in  which  the  fuel  is  embedded  is  infinite  spatially.   Practical 
experimental  assemblies  for  neutron  wave  propagation  experiments,  however, 
tend  to  be  quite  small  in  comparison.   Consequently,  it  was  necessary  to 
start  from  the  basic  equations  including  the  finiteness  of  the  feasible 
experimental  assemblies. 


CHAPTER  I 


USE  OF  A  ONE -GROUP  ONE -DIMENSIONAL  GREEN'S 
FUNCTION  IN  HETEROGENEOUS  MEDIA^ 


The  problem  to  te  considered  is  the  calciilation  of  the  transmission 
of  a  thermal  neutron  wave,  according  to  one-group  diffusion  theory,  in  a 
simple  heterogeneous  geometry.   The  geometry  chosen  is  such  that  it  permits 
a  practical  experimental  setup  and  also  results  in  an  analytical  form  reduc- 
ible virtually  to  a  one -dimensional  problem. 

Figure  1  shows  the  geometry  of  the  problem.   The  semi-infinite 
parallelepiped  of  moderating  material  extends  from  x  =  -a  to  x  =  a,  y  =  -b 
to  y  =  b,  and  z  =  0  to  z  =  oo.  The  transverse  dimensions  a  and  b  are  assumed 
to  include  the  diffusion  theory  extrapolation  distance.   The  perturbation 
in  the  medium  is  taken  to  be  in  the  form  of  very  thin  sheets,  the  kth  sheet 
extending  from  x  =  -x^  to  x  =  x^  and  y  =  -b  to  y  =  b  and  of  thickness  h, 
centered  aro'und  z  =  zj^.   The  system  is  driven  by  a  sinusoidally  modulated 
soirrce  of  thermal  neutrons  placed  at  the  face  z  =  0. 

The  time -dependent,  one-group  diffusion  equation  describing  this 
system  is 

i  ^^=  {iP^   -      [I^.    BZjr)]}  .(r,t)   ,  (1.1) 


"""The  solution  obtained  in  this  chapter  is  an  extension  of  work 
originally  performed  by  R.  S.  Booth. 


Fig.  1.   Geometry  for  the  One-Group,  One-Dimension  Green's  Function 
Calculation. 
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the  neutron  speed, 

the  spatial  and  time -dependent  neutron  flux, 

a  position  vector, 

time, 

the  diffusion  coefficient  (units  of  length), 

the  Laplacian  operator, 

the  macroscopic  cross  section  of  the  moderating  media, 

the  space -dependent  macroscopic  cross  section  of  the 

absorbing  sheets. 


The  space -dependent  perturbation  may  be  ^.■nritten  as 


M 
6zjr)  =  ^  6Z^^ 

10=1 


u(x  -f  x^  -  u(x  -  Xq) 


u  I  z  -  Zj^  +  I 


u  z 


h 


(1.2) 


where  6Z    =  the  macroscopic  cross  section  of  the  kth  foil,  u(y)  =  the 
a,  X  — 

unit  step  fxmction,  and  M  =  the  total  number  of  foils.   If  h  is  taken  to  be 
s\ifficiently  small,  the  first  two  terms  of  the  Taylor's  series  expansion 
for  the  unit  step  function  in  z  may  be  used  in  Eq.  (1.2),  reducing  it  to 


6Z  (r)  =  h 
a  — 


u(x  +  Xq)  -  u(x  -  Xq) 


M 


5Z  -  5(z  -  z,  )   , 
a,k  ^     k   ' 


(1.5) 


k=l 


where  5(y)  is  the  Dirac  delta  function. 


It  shoiald  be  noted  that  if  the  sheets  contain  fissionable  material^ 

5Z    is  replaced  by  62  ,  -  v52  ,,    where  v  is  the  neutron  multiplicity 
a  •  X  a^  K.  ? 

and  52  ,  is  the  macroscopic  fission  cross  section  of  the  kth  foil. 

I^K.  — 

The  sinusoidally  modulated  neutron  source  is  assumed  to  be  isotropic 
with  a  Maxwellian  energy  distribution  and  can  be  represented  by 


S(r,t)  =  S  +  Re 


S(r) j    e 
z-o 


itot 


(1.^) 


where  i  =  v-1^  '^  i^  "the  angular  frequency  of  the  modulation  and  S  is  the 
steady  component  of  the  source  present  at  z  =  0.   If  the  neutron  flux  is 
assumed  to  be  separable  in  space  and  time^  the  oscillating  part  of  the  flux 
can  be  represented  as  (l) 


(t>(r^t)  =  Re 


(J)(r) 


iwt 


(1.5) 


Equation   (l.5)    may  be   introduced  into  Eq.    (l.l),    resulting  in 


^    - 


2+52   (r) 
a a  — 

D 


ICO 


*(£)  =  7PR  *(£)     ^ 


vD 


(1.6) 


or^  in  slightly  different  notation^ 


^ 


.(r)  =  5L(x,z)  (t)(r)  , 


(1.7) 


where 


a„  =  v2 


a 


D^  =   vD 


6L(x,z)  = 


52  (r)/D 


(1.8) 

(1.9) 

(1.10) 
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The  boundary  conditions  will  be  taken  to  be 


>(r) 


=   0(r) 


x=+a 


y=+b 


lim      .  /    V        ^ 
<t>(rj   =   0 
z-*oo        — 


(1.11) 


In  addition,  at  the  source  plane, 


d<t)(r.t) 
dz 


z=o 


=  I  Re 


S(r) I     e 
z=o 


iwt 


(1.12) 


Using  Eq.  (I.5),  this  reduces  to 


D 


dt)(r) 


az 


=   |s(r)| 
z=o         z=o 


(1.13) 


The  spatial  flux,  'J'(r),  may  be  represented  by  an  expansion  in  the 
complete  set  of  transverse  eigenfuactions  defined  by 


(V^  +  B^,  )  <t>,  (x,y)  =  0 


il.lh) 


and  by  the  boundary  conditions  at  x  =  +a  and  y  =  +b  plus  the  conditions  of 
symmetry  around  x  =  0  and  y  =  0.   Normalized,  these  eigenf unctions  and  their 
related  eigenvalues  are 


1      (2i  -  l)«x     (2m  -  l)7ty 

-fl  \--5j/  -  COS  -^ — - —  cos  -i — — '—^ 

M^    ''^'         rr-        2a  2b 


n^(^^y) 


(1.15) 


and 


jLjfim 


(2i  -  l)ji 
2a 


-,2 


(2m  -  l)it 
2b 


(1.16) 


where  Z   and  m  may  take  on  all  integral  values, 


The  expansion  of  0(r)    is  given  by 


(1.17) 


Substituting  Eq.  (I.I7)  into  Eq.  (I.7)  gives 


^ 


a  +  iw" 
o 

57 


£,Ta 


a,       u 
Operating  on  Eq.  (I.I8)  with   /   dx   /   dy  *   (x^y)  results  in 


-a 


(1.19) 


where 


p2  =b2 


Gq  +  ioj 


pq    J-pq     D^ 


(1.20) 


and 


a     b 


^^pqM^^)  ^    J    ^"^    J    ^^   Vq^^^y^  ^^^"^")  *-^"^y) 


M' 


-a    -b 


(1.21) 


Before  continuing  further  with  Eq.  (I.I9),  note  that  the  two  boundary 
conditions  in  z  have  yet  to  be  applied.   In  order  to  make  use  of  the  source 
boundary  condition^  S(r)  must  be  expanded  in  the  same  transverse  eigenfunc- 
tions.   That  is^ 


S(r) 


10 


where 

a     b 
S 


dx   /dyS(r)|    ^Jx^y)   •  (1-25) 

^  '7  — r. 


-a    -b 


It  should  be  pointed  out  that  this  expansion,  as  well  as  the  flux  expansion, 
assumes  symmetry  aroimd  the  x  and  y  axes.   If  one  uses  a  source  containing 
asymmetric  components,  *   (x, y)  must  also  include  sine  terms.   Inserting  the 
expansions,  Eqs.  (I.17)  and  (1.22),  into  the  source  boimdary  condition, 
Eq.  (1.13),  gives 


d<D.  (z) 


D   '-^ 


dz 


S 


^     .  (1.2U) 


z^o 


The  boundary  condition  at  z  =  00  becomes 


^^"^  n  (z)  =  0  .  (1.25) 

z^oo   ^m 


Now  returning  to  the  reduction  of  Eq.  (I.19),  substitute  Eqs.  (I.5) 

and  (1.10)  into  Eq.  (l.2l).   The  y  integration  may  be  performed  immediately, 

yielding  6  ,  the  Kronecker  delta.  Also  performing  the  x  integration  gives 

the  result 

M 

5L  ,  (z)  =  ^  5   5F  ,  )   5Z  -  6(z  -  zj   ,  (l.26) 

pqim      D  qm   p^  /_,    a,k        k   ■* 

k=l 


where  5F  .  is  given  by 
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(i  -  p)itx^        (-2  +  P  -  1)«Xq 
sin  sin 

^^Pi  =       {z   -  p)n    -^  — TTTT^TTT—  '  ^  ^  P 


(1.27) 
(2i  -  l)nx^ 

--^    (2i-l)n     .   ^  =  P   • 

In  order  to  simplify  Eq.  (1.27)  fiirther,  let  us  assume  that  x^  is  suffi- 
ciently less  than  a  that  the  first  term  of  the  Maclaurin's  series  may  be 
used  for  the  sine  terms.   Equation  (1.27)  then  becomes 

2x 

6F ,      for  all  Z   and  p.  (1-28) 

pi  -  a 

Wow  inserting  Eq.  (1.26)  and  Eq.  (1.28)  into  Eq.  (1.19)  gives 


Ldz2     PI- 


M 
2hx 


<t>   (z)  = )  )   52   ,  5(z  -  z  )  <{,   (z)   .         (1.29) 

pq^      a    /_,   Zl^    a,k        k   iq 


H,        k^l 


It  may  be  noted  from  the  above  equation  that  the  one -dimensional  problem 
turns  out  to  be  not  really  one  dimensional^  in  the  sense  that  coupling  re- 
mains between  the  spatial  harmonics  associated  with  the  x  direction. 

In  order  to  solve  Eq.  (1.29),  it  is  convenient  to  convert  it  into  an 
integral  equation.   Substitute  %    for  z  and  rewrite  as 


^     (I)  =  F  (I)  ,  (1.50) 

v^  pq 


where  „ 

2hx 


F  ii)^     ^    )       )      6Z  ,5(1  -  z,  )  <t),  (I)   .  (1.51) 

pq^^'^    Da   ^  /_j   a,k   "=    k   iq^' 

£     k=l 
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Now  define  the  Green's  function  G  (z,^),    satisfying 


a^      2 
P 

2   pq 


sr 


v^^^)  - 


D 


(1.32) 


and  the  homogeneous  forms  of  the  boundary  conditions  satisfied  by  <t)   (z) 
Multiply  Eq.  (l.JO)  by  G   (z,t)  and  Eq.  (1.32)  by  ^  (i) .      Subtract  the 

00 

two  results  and  operate  with   /   d|  to  obtain 

o 


00 


d^  iG(z,0 


d^»(0 


oo 


-  j   d^  F(|)  G(z,|)  = 


»(z) 
D 


(1.33) 


The  subscripts  have  been  dropped  temporarily  for  economy  of  notation. 
Integrating  the  first  integral  by  parts  results  in 


*(.)  =  b{g(.4)^-*(5)^^} 


oo 


oo 


^  J  ^^    al       dl 


00 


+  D     d| 


dO(|)   SG(z,|) 


oo 


D  /   d|  F(5)  G(z,a; 
o 


(1.3^) 


The  first  two  integrals  cancel,  both  *  and  G  go  to  zero  as  |  ->  oo,  and  the 
gradient  of  G  is  zero  at  ^  =  0.   Therefore,  the  desired  integral  equation 


IS 


00 


$(z)  =  -  DG(z,|)  ^ 


-  D    J      d^   F(5)  G(z,5) 
1=0      o 


(1-35) 
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The  problem  now  is  to  determine  the  Green's  function  defined  by 
Eq.  (1.32).   The  solution  to  the  homogeneous  form  of  Eq.  (I.52)  is 


pz   „  pz 


G(|,z)  =  Ae  ^"  +  Ce^"   , 


(1.36) 


where  |  and  z  have  been  switched.   The  particular  integral  of  Eq.  (I.52)  may 
be  found  by  the  technique  of  the  variation  of  parameters.   Let  the  constants 
in  Eq.  (I.36)  depend. on  z   (they  also  depend  on  i,    of  course,  but  this  may 
be  ignored  at  this  point) : 


pz 


G(|,z)  =  A(z)  e"^"  +  C(z)  e^"  . 


(1.57) 


Differentiation  of  Eq.  (I.37)  gives 


^G(g,z) 
Sz   ^ 


.  -pz     „  pz   dA  -pz   dC  pz 
-  p  Ae    +  p  Ce*^  +  --  e  "^  +  -—  e"^ 
^         ^        dz        dz 


(1.38) 


Choose  A  and  C  such  that 


dA  -pz   dC  pz    _ 
dz        dz 


(1-39) 


Then  differentiate  again: 


az2 


dA  -pz     dC   pz 
p-—  e"^  +p— -e 
dz  dz 


(1.^0) 


Inserting  Eqs.  (l.^O)  and  (1-37)  into  the  original  differential  equation 
for  the  Green's  function  results  in 


A 


2  -pz 
P  e 


2  -pz 
P  e 


+  C 


2  pz    2  pz 
p  e   -  p  e 


+ 


dA 
dz 


pe 


-pz 


dC 
dz 


pe 


,PZ 


S(z  -  I) 
D      ' 


(l.i+1) 


ll^ 


or 


-pz  dA     pz  dC     S(z  -  f)  1^    ,  ^v 


Equations  (I.59)  a-^i'i  (l-^2)  are  two  eqiiations  in  dA/dz  and  dC/dz  which  may 
be  solved  simultaneously  to  give 


dA  _  5(z_:_|l_e^  (IW) 

dz  -    2pD  u-^^; 


and 

dC_  _  5(z  -  g)  e-P^  ,   j^^. 

dz  "        2pD       •  ^^-^^^ 

Integrating  Eqs.  (lAj)  and  (l.Uij-)  gives 

A  =  ePV2pD  +  B  (1.^5) 

and 

C  =  -  e-^VspD  +  E  ,  (1-^6) 

where  B  and  E  are  constants  to  be  determined  from  the  boundary  conditions. 
Using  Eqs.  (l.ij-5)  and  (I.U6)  in  Eq.  (I.56)  gives 

-DZ     oz   eP^^"^^    „-p(|-z) 
G(g,z)  .  Be  P^  .  EeP^  H-  ^-^-^  -   ^^^^—  (l.^T) 


or 

G(g,z)  .  Be-P^  +  EeP"  +  ^^^  ^^^  "  "^   .  (1A8) 

pD 


The  one -dimensional  Green's  fimction  has  a  discontinuity  in  slope  at  the 
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source  pointy  z  =  | .  Actually^  this  does  not  need  to  be  known  "beforehand, 
since  straightforward  application  of  the  only  two  boundary  conditions  now 
available  leads  to  the  solution  G  =  0.   Then,  in  order  to  obtain  a  solution, 
before  applying  the  boundary  conditions,  Eq.  (l.^8)  must  be  divided  into 
two  parts,  one  which  will  apply  for  z  <  |  and  the  other  for  z  >  | .   Now, 
two  more  boimdary  conditions  are  needed,  which  may  be  found  by  integrating 
Eq.  (1.32)  over  a  small  region  including  the  source  point: 


;+e 


dz 


|-€ 


LdZ-^ 


;+€ 


G(|,z)  = 


dz 


5(z  -  g) 


D 


(l.i+9) 


l-e 


Assuming  that  the  second  derivative  is  more  singular  than  the  function  it- 
self, one  obtains,  in  the  limit  of  small  e, 


Sz 


§z 


1 
D 


(1.50) 


l+e 


|-€ 


Therefore  the  two  boimdary  conditions  are  continuity  of  G  at  z  =  |  and 
Eq.  (1.50),  which  gives  the  magnitude  of  the  discontinuity  in  the  first 
derivative. 

The  two  parts  of  the  solution  will  be  denoted  as 


L  -pz    L  pz    sinh  pjj    -   z) 


g"(|,z)  =  B"  e"^"  +  E"  e^^  + 


pD 


.   z  <  ^   , 


(1.51) 


and 


G^|,z) 


R  -pz    R  pz   sinh  p(g  -  z) 


B  e 


+  E  e 


pD 


z  >  I 


(1.52) 


The  z  =  0  boundary  condition  is 


dG^^,z) 


dz 
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PB^  +   pE^   -  ^^^  =    0 


z=o 


At  large  z,    the  boundary  condition  is 


(1.53) 


lim  _R/.   s    lim 
G  (|,z)  = 


'r     pz   e-P^  e  P" 


=  0 


(1.5^) 


R 


Equation  (1-5^)  immediately  determines  E  : 


E^=^ 


-Pi 
2pD~ 


(1.55) 


The  contin\iity  condition,  at  z  =  |,  leads  to 


B 


h-^^   +  E^eP^  =  B^e-P^  +  ^\'' 


(1.56) 


Equation  (I.50)  becomes 


pB^-P^  .   PE^P^  -  i  +  pB^e-P^  -  pE^eP^  .  i 


1 
D 


(1.57) 


Using  Eq.  (1-55),  Eqs.  (1-55),  (1-56),  ajid  (l.5T)  may  be  put  in  matrix  form 


as 


-  P 
-Pi 


P 
,Pl 


B 


cosh  p| 
D 


-Pi 


-pf      Pl      -Pl 
pe  ^^    pe*^^    pe 


E 


2pD 


B 


R 


\    1- 

\       2D 


(1.58) 


The  three  solutions  to  Eq.  (I.58)  are 


IT 


L_   sinh  Pi  (1.59) 


■Pi 
E  = 


pD   ^ 


-Pi 
B  = 


R   e  -  (1,61) 


2pD 

Upon  substituting  Eqs.  (1-55);  (1.59),  (I.60),  and  (I.61)  into  Eqs.  (I.51) 
and  (1.52)  and  rearranging^  one  obtains 

G^|,z)  .£!L£o^lL^,   .<!   ,  (1.62) 

and 

G^|,z)  =£^1|°S1L^,   z  >|   .  (1.65) 

It  may  be  noticed  that  Eqs.  (1,62)  and  (I.65)  define  a  function 
which  is  symmetric  with  respect  to  interchange  of  the  coordinates  of  the 
source  and  observer  locations  (z  and  i) ,    thus  satisfying  one  of  the  nec- 
essary criteria  for  a  Green's  function  (9).   Equation  (I.63)  may  now  be 
substituted  into  Eq.  (1-35),  giving 


(t)(z)  =  -  D 


''      -  -D  f€^^^^L^Hi)    dl 


PD   d^ 


1=0     o 

00 


.  D  /■  e  P^  cosh  ^  ^(^^  ^^   _  (1,6,,) 

J       pD 


z 


Reinserting  subscripts  and  using  Eq.  (1.2^)^  one  obtains 
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-p  z 

o    pq.    -p  z   z 

S   e  -^^       pq    n 

<D   (z)  =  -^ r^ - /  d|  cosh  p  i   F     (t) 

•^^  pq         pq    ^  jr-i   JT-TL 


cosh  p  z  p      -P   I 

pq   /   d|  e  P^  F  (0   .       (1.65) 

pq    '^ 


Since  F  (0  contains  5  functions  in  |,  Eq.  (l.Jl)  may  be  inserted  into 
pq 

Eq.  (1.65)  and  the  integrals  in  |  performed^  leading  to 

S   e'^'P'l^   2hx^    r-p  z  ^  f 

0   (z)  =  ^ea-^ ^  Je  P'l   )   )  6Z  ^  cosh  p   z^  *   (z,  ) 

pq^  ^      2p  D     p  Ea   t       Z_j  Z^    a-^k      "^pq  k  iq'  k' 
pq.      P<1    ^       ^  l^T^ 

M 

.  cosh  p^^z  1       y^  5Z^^^  e"'P-l'>^  *^^(z^)}  ,    (1.66) 

£   k=M'+l 


where  M'  is  determined  by  z  ,  <  z  <  z  ,   .   If  z  equals  one  of  the  z  's^ 

it  can  be  taken  to  be  either  z,,,  or  z,,,  ^    since  the  Green's  function  is 

M'     M'+l 

continuous  at  the  source  point. 

In  order  to  obtain  a  solution  to  this  system^  the  original  expansion 
of  •t'Cr),  Eq.  (1.17)^  must  be  truncated.  Assume  that  the  desired  approxi- 
mation is 


*(£)  -  /_.  A  ^q^^^  ^iq^""^^^   •  ^^'^^^ 

£^1   q=l 

For  each  value  of  q,  a  set  of  LM  simultaneous  equations  may  be  obtained 

from  Eq.  (I.66)  by  writing  it  for  each  of  the  L  modes  and  for  z  =  z  . 

z„,  ...  Z|«.   This  set  determines  each  of  the  LM  values  of  0.  (z^  ). 
2-'      i--i  iq^  k 
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Writing  these  values  as  a  vector^  where  the  superscript  T  indicates  the 
transpose^ 


$^  = 


W\^'    %(^2)^ 


^q^^M^^  %(^1^^  •••  %^V 


(1.68) 


the  set  of  equations  may  be  put  into  matrix  form.   That  is^ 


A 


1   L 


$  =  B 

q  1 


(1.69) 


where  [l]  is  the  identity  matrix,  the  vector  B  contains  the  constant 

q. 

(source)  terms,  and  the  matrix  [A  ]  contains  the  Green's  functions.   B 

q  q 

may  be  written  as 


q. 


a   e 

.  1^ 


-p  z 
^Iq  1 


a,  e 

1<1 


''^iq's 


a^  e 
Lq 


-D   Z 


(1.70) 


where 


a 


£q 


iq   2p^^D 


(1.71) 


The  [A  ]  matrix  may  be  partioned  into  L^  submatrices  which  are  each  M  by 
M.  The  n^p  element  of  the  I^J  submatrix  is  given  by 


(A^  t-)    =  P^  (z  )  cosh  p^ze   '^'^.n>p, 
I, J  n,p    Iq  p'       Iq  P  '         ' 


(1.72) 


or 


-P-r  Z 

(a   )    =  P^  (z  )  cosh  p^z  e       ,    n<p   , 
I,J  n,p   "^Iq'  p'       "^Iq  n         '  ^      ' 


(1.73) 


where 
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2hx  52 
^Iq^^p^     Pj^  Da 


(1.7^) 


Solving  Eq.  (I.69)  Q  times  by  inverting  the  matrix  [A  ]  +  [I]  gives 


the  $  's: 

q. 


[A  ]  +  [I] 

q. 


(1.75) 


which  may  be  inserted  into  the  s-ummations  in  Eq.  (I.66).   The  results  of 
Eq.  (1.66)  are  then  inserted  into  Eq.  (I.67)  to  give  the  desired  approxima- 
tion to  the  spatial  flux  <l>(r). 


CHAPTER  II 


USE  OF  A  ONE-GROUP  TWO-DIMENSIONAL  GREEN'S 
FUNCTION  IN  HETEROGENEOUS  MEDIA 


The  problem  of  the  previous  chapter  will  now  be  extended  to  two 
dimensions  by  allowing  the  perturbations  to  be  placed  anywhere  in  the  x-z 
plane  contained  in  the  moderating  material.  Equation  (l.l)  will  still  apply, 
but  now  5Z  (r)  will  be  taken  to  be  a  general  function  of  x  and  z.   Eqimtion 
(1.7),  with  its  bo-undary  conditions  given  in  Eqs.  (l.ll)  and  (I.I5),  will 
be  taken  as  the  starting  point,  leaving  6L(x, z)  general. 

In  this  case  the  spatial  flux  <t)(r)  will  be  expanded  in  the  single 
set  of  transverse  eigenfunctions,  defined  by 


m 


dy^ 


%(y)  =  0  (2.1) 


and  the  boundary  conditions  0  (+b)  =  0  and  'i>    (y)  =  *  ( -y)  •   The  eigenfunc- 
J  m  —  mm 

tions  are 


V,,  -   ^  cos  (5L^^  ,  (..2) 


and  the  eigenvalues  are 

(2.3) 
Using  the  following  expansions  of  the  flux  and  the 
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■r2 

(2m  -   l)n 

£i        — 

m 

2b 

where  m  =   1,2,3,  .  . 

source 

9 
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*(r)  =  )   V^^^"*  V^^'' 


m 


and 


S(r)  .  ^  S^(x)  .^(y)   , 


m 


(2.U) 


(2.5) 


in  Eqs.  (l-T)  and  (I.I5)  gives 
a  +  iw 


^ 


D 


^  V^^^^)  V(y)  =  6L(x,z)  ^  V(x.^)  \(y)    .  (2.6) 


m 


m 


and 


-  D  ^  y  0  (x,z)  My)     =  i  y  S  (x)  <t>  (y) 


m 


z=o      m 


(2.7) 


u 
Operating  with   /  dy  <J>  (y)  gives 


52      a^ 


ox^   dz^ 


D^(x,z)  =  5L(x,z)  -f^Cx^z)   , 


(2.8) 


where 


a  +  iw 


p2  ,  B^  +  ^ 
n    n     D 


(2.9) 


and  the  corresponding  boundary  condition 


-  D 


a*  (x.z) 
n  ■' 


=  i  S„(x) 


z=o 


(2.10) 


25 
The  remaining  boimdary  conditions  given  in  Eq.  (l.ll)  still  apply. 

It  is  evident  that  this  system  may  be  treated  similarly  to  the  prev- 
ious problem  if  a  two-dimensional  Green's  function  can  be  found  which  is 
the  solution  of 


af__^  af_ 


dx^     ^2 


n 


G^^^^^^^o^^o^ 


6(x  -  Xq)  5(z  -  Zq)  , 


(2.11) 


with  boundary  conditions 


SG 


G  (.+a,z;x  ,z  J  =  ^; — 
n^-  -^   o^  o'   dz 


^^"^  G  =  0 
z-»co   n 


z=o 


(2.12) 


This  Green's  function  may  be  found  by  expanding  in  the  ei gen functions  of 


dx^ 


k2 


\l;^(x)  =  0;   tj^(+a)  =  0 


(2.15) 


These  are 


,  /■  \    .               (2i  -  l)nx  ,  ^    .   iitx 
^^U)    =  A_g  cos  -^ +  B^  sm  — 


(2.11^) 


Note  that  k  assumes  two  different  values  for  the  two  different  eigenf unc- 
tions^ the  sine  and  the  cosine.   The  expansion  of  G  (x, z;x^^z^)  is 


G^(x,z;x^,z^)  =  ^  C_g(z,z^)  k^U^)    cos  (2-g  -^l)  nx  ^^^^^^^  ^.^ 


ijtx 


(2.15) 


Inserting  Eq.  (2.15)  into  Eq.  (2.11)  gives^  with  arguments  dropped  for  con- 


venience. 


21+ 


2  2 

„      f{2£    -   l)n\  {2£    -   l)«x       „     j£n\        .      £nx 


i 


(2i    -  l)itx       „        .      inx 


p|  /^. 


[2Z    -    l)rtX        ^         .      i:rtx 
h   ''°'  2i ^  ^i    ^"^  — 


-  5(x  -  Xq)  5(z  -  z^)   .   (2.16) 


a 


Operating  on  Eq.  (2.l6),  in  turn,  with   /   cos  -^ — ^—^ — - —  dx  and  with 


a 


-a 


sm  '^ —  ox  gives 


-a 


C.  A. 


|(2j^LA)ill  ^  ^  _J.  A.a  -  p2  C.A.a  =  -  6(z  -  zj  cos  ^^J  "  ^^  ^""^   , 
I   2a    J      ^^2   J     ^  J  J       ^     o^         2a      ^ 

(2.17) 


and 


r.  ^2   a^c. 

-  C.B.  i^^  a  +  ^  B.a 

J  J  la  J     ^^2   J 


p2  C.B. a  =  -  6(z  -  z^)  sin  

n  J  J  O''      a 


(2.18) 


These  two  expressions  may  be  rewritten  so  as  to  have  the  left-hand  sides 
fionctions  of  z  and  z  only  and  the  right-hand  sides  functions  of  x  only. 
They  become 

(2j  -  l)TtX^ 


3z2   L  n   L   2o 


C.    cos — 

J  2a 


-  6(z  -  Zq) 


(2.19) 


a  A. 
J 


and 
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a^c 


dz2 


_1  _ 


p2  H^ 

n    a 


Jrtx^ 


sm  — 


a 


-  &(z  -  Zn) 


(2.20) 


a  B. 


Obviously  the  C.'s  in  the  above  two  equations  are  not  the  same;  so  the  C. 
J  J 

will  be  replaced  by  Al  in  Eq.  (2.19)  and  by  B!  in  Eq.  (2.20).   Notice  also 

3  J 

that  both  sides  of  Eqs.  (2.19)  and  (2.20)  are  equal  to  constants  and  no 
loss  in  generality  occurs  if  this  constant  is  taken  to  be  unity  in  each 


equation.   With  this  step,  A.  and  B.  may  be  obtained  immediately: 


A.(x^)  =  icos 


(2j  -  l)nx^ 
2a 


(2.21) 


and 


B.(xo)  = 


1   .   J"^o 

—  sm  

a      a 


(2.22) 


and  with  the  definitions 


(2.23) 


^nd 


a 


(2.2U) 


the  two  equations  for  A I  and  Bl  are 


L^z^ 


-  M- 


A'.(z.z  ) 


-   S(z  -  Zq) 


(2.25) 


and 
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a^ 


ri^ 


B1(z,Zq)  =  -  5(z  -  Zq)  .  (2.26) 


Laz^ 

The  related  boundary  conditions  are 


SB! I      ..        .. 
^  ^     ^  lim  a:  =  ^^"^  Bl  =  0  ,  (2.27) 

dZ  I        Z-KX3   J    z-*oo   J        ^ 


z=o       z=o 


plus  the  function  continuities  and  first-derivative  discontinuities  at 
z  =  Zq.   These  functions  are  virtually  identical  to  the  one-dimensional 
Green's  function  developed  in  Chapter  I  and  can  therefore  be  immediately 
written  down,  being 

„  ,  /     V    e     cosh  LiZ        , 
A  (z^Zq)  =  ^  ,    ^  <^o 


e    cosh  LIZ 
^  o 


z  >  z^   ,  (2.28) 


>         "    -    "o 


and 


_,/    V   e     cosh  riz 


e  '  cosh  Tiz 
'  o 


,   z  >  Zq   .  (2.29) 


Collecting  Eqs.  (2. 21),  (2.22),  (2.28),  and  (2.29)  into  Eq.  (2.15),  the 
complete  Green's  function  is 


27 


^n^^'^^^^'o'^o)  =^ 


e cosh  \iz 

L      M-a 


{2£    -   l)rtx 


cos 


2a 


o     (2i  -  l)nx 
"°^  2^ 


^ 


e     cosh  nz   .     o   .   ijtx 

+  *—  sm  sm 

Tja  a        a 


>        2  <  ^o 


(2.30) 


and 


G  (x,z;x^.z^)  = 
n  '   o'  o 


I 


^'       ^°"^  ^^o     (2i-l)nx^      (2i  -  l)nx 


|ia 


cos 


2a 


cos 


2a 


e  '  cosh  t^Zq     -^^^o 


+ 


T^a 


sm 


a 


sm 


a 


z  >  z 


(2.51) 


Similarly  to  the  procediore  followed  previously^  this  Green's  function 

may  be  used  to  formulate  an  integral  equation  for  <t)  (x,  z).   x  and  z   are 

interchanged  with  x  and  z  in  Eqs .  (2.8)  and  (2.11).   The  difference  between 

G  (x^z  ;x^z)  times  Eq.  (2.8)  and  $  (x  ^z  )  times  (2.11)  is  taken.   The 

GO      a 
result  is  operated  on  with   /   dz^   /  dx^  to  give,  with  arg-uments  and  inte- 


-a 


gration  limits  dropped. 


^2* 


B2* 


'^^o   /    ^o^n 


n 


< 


+        dz_        dx     G 


n 


S^G 


dz      /    dx„    4) 
o,/        o     n 


ax§ 


a^G 


dz^    /    dx^   * 


n 


dz      /    dx„  G     5L  t-      +   <t>    (x.z)       ,  (2.52) 

o,/        On  n         n-" 


In  the  terms  containing  partials  with  respect  to  x,  an  integration  by  parts 
in  X  may  be  done,  and  similarly  with  the  z  partials.   All  remaining  surface 
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integrals  on  the  left-hand  side  of  Eq.  (2.52)  cancel  and  most  of  the  line 
integrals  are  eliminated  "by  the  boundary  conditions.  The  remaining  terms 
give 

0^(x,z)  =  -  /  dXg  G^(Xq,Zq;x,z)| 


5z^ 

a  zcfo  Zo=o 


00      a 

J     ^^o    J   ^0  ^n^^'o^^o''''^)  S^(^o^^o)  V^^o^^o^   •         (2.53) 


-a 


Equations  (2.10)^  (2.5O),  and  (2.51)  are  to  he  substituted  into  Eq.  (2.35), 
along  with  the  desired  form  of  5L(xo^2,o),  in  order  to  obtain  the  integral 
equation  for  $  (x, z). 

If  the  form  of  6L(x, z)  of  Chapter  I  is  used  in  Eq.  (2.55)^  i"b  is 
straightforward  although  rather  tedious  to  show  that  Eq.  (2.33)  reduces  to 
Eq.  (1.66).   To  do  this,  the  following  expansions  must  be  used: 

<t>  (x.z)  =  )   (D.  (z)  <D.(x)  (2.3U) 


and 

S  (x)  =  >   S.   4).(x)   ,  (2.35) 


where 

♦  .UJ.^eosM^jiis  .  (2.56, 

V  3- 
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These  expaxisions  are  justified,  since  if  6L  is  symmetric  in  x  [and  it  is 

assumed  as  before  that  S(r)  is  symmetric  around  the  x  and  y  axes],  then 

<t>   will  also  be  symmetric  in  x.   It  is  then  a  matter  of  performing  the  Xn 
n 

a 
and  Zq  integrations  in  Eq.  (2.55)  and  operating  with   /  ^—  cos     „ dx 


to  obtain  Eq.  (I.66) 


CHAPTER  III 


USE  OF  M   AGE -DIFFUSION  TWO-DIMENSIONAL  GREEN'S 
FUNCTION  IN  HETEROGENEOUS  MEDIA 


The  development  of  this  chapter  is  similar  to  that  of  the  preceding 
chapter  except  that  the  theory  is  to  be  the  continuous  slowlng-down  model 
with  a  thermal  group.   The  geometry  is  the  same  semi-infinite  moderating 
medium,  with  the  perturbations  now  being  small  fuel  rods  which,  as  before, 
run  the  full  length  of  the  y  axis.   This  model  is  similar  to  the  Feinberg- 
Galanin  treatment  of  heterogeneous  reactors  (7)  and  (8),  with  two  major 
exceptions.   In  this  study  time  dependence  is  included  and  the  finiteness 
of  the  medium  is  taken  into  account.   This  will  result  in  a  finite  medium 
Fermi -age  kernel,  or  Green's  function,  as  well  as  the  corresponding  diffu- 
sion kernel  obtained  in  the  previous  developments  of  this  work. 

For  lethargies  below  thermal,  the  applicable  equation  is 
-  D(u)  ^   1>(r,u,t)  +  Z  (u)  <t(r,u,t) 

1    d»(r,u,t)  bq{r,n,t) 

=   -  —T—\ v;: +  S(r,u,  t) T ,  15- 1; 

v(u)     ot         —'    '  ou       ' 

where  the  nomenclature  is  similar  to  that  previously  used  except  that  now 

all  qimntities  are  lethargy-dependent  and  S  and  q  represent  the  neutron 

source  and  slowing-down  density,  respectively.   It  has  been  assumed  that 

Z  (u)  for  u  <  u  ,  where  u  is  the  thermal  neutron  lethargy,  is  independent 
a  s'        s 

of  position.   This  is  valid  as  long  as  the  fuel  rods  are  small. 
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The  relation  between  the  flux  and  slowing-down  density  applicable 
to  an  infinite  medium  will  be  assumed^  that  is^  that 

q(r,u,t)  =  i    Z^(u)  *(r,u,t)   ,  (5-2) 

where  |  -  the  average  logarithmic  energy  decrement^  and  2  (u)  =  the  total 
macroscopic  cross  section  in  the  moderating  media. 

The  thermal  group  is  described  as  before^  except  for  the  addition 
of  the  source  of  thermal  neutrons  from  the  slowing-down  density.   This 
relation  is 


D  ^  +  2   +  —  It 

s       as   V   ot 
s 


^'Jr.t)  =  q(r,u  t)  -  A(r,t)   ,  (3.5) 


where  A(r, t)  represents  the  absorption  in  the  fuel  rods  and  the  subscript 
s  indicates  that  the  quantity  is  evaluated  at  the  thermal  neutron  energy. 

In  order  to  define  the  quantities  S(r, u^t)  and  A(r^t),  it  will  be 
assumed  that  (l)  all  fission  neutrons  are  bom  at  lethargy  u  =  0^  (2)  all 
fuel  rods  are  lines^  the  kth  rod  being  representable  as  S(x  -  x^)  S(z  -  z  ), 
and  (3)  the  driving  source  is  located  at  z  -  0,  producing  only  thermal  neu- 
trons and  is  of  the  form  Re[S(x,y)  6(z)  e    ].   As  before,  since  the  source 
is  present  only  at  the  boundary  z  =  0,  it  will  be  included  as  a  boundary 
condition  rather  than  being  included  in  Eq.  (5-5) • 

Since  this  treatment  is  an  extension  of  the  Feinberg-Galanin  treat- 
ment, some  of  Feinberg's  notation  will  be  adopted.   Define  the  thermal 

constant  7,  as  the  ratio  of  the  total  net  current  of  thermal  neutrons  into 
k 

the  kth  rod  to  the  value  of  the  thermal  neutron  flux  at  the  surface  of  the 

kth  rod.   0  (r,  ,t)  will  be  used  for  the  thermal  neutron  flux  at  the 
—         s  — k^ 
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surface  of  the  kth  rod.   The  neutron  yield  per  thermal  neutron  absorption, 
T],   will  "be  ass-umed  to  be  the  same  for  all  rods.   With  the  above  assumptions 
and  definitions,  S(r, u, t)  and  A(r,t)  become 

S(r,u,t)  =     ^  6(u)  ^-^  T]  7^   %(£j,.*)  5(£  "  £1,)  (5-^^) 

1  ^ 


and 

A(r,t)  =    /_^    y^   *s^-k'*^  ^^-  "  -k^   •  ^^"^^ 

k 

The  ratio  of  the  scattering  to  total  cross  sections  at  the  source  lethargy 
in  Eq.  (5.^)  represents  the  probability  of  a  neutron  starting  the  slowing- 
down  process. 

Due  to  the  separability  of  space  and  time  for  the  flux,  slowing- 
down  density,  and  source  term,  Eqs.  (5.I),  (3.2),  and  (j.U)  may  be  combined 
to  give 


-  D(u)  ^  +  Z  (u)  +  4\  +  I  Z^(u)  I- 
a      v(u)      t    du 


"1 

q(r,u) 


ZjO)    ^ 

=  ^  ^^"^  rroy  ^  I  ^^-^  ^k  ^s^^k^  ^^^  -  ^k^   •  ^5.6) 


Similarly,  Eqs.  (5- 3)  and  (3-5)  become 


D  ^  +  2   +^ 
s       as   V 

s 


t)Jr)  =  q(r,Ug)  -  ^^    y^  ^^{r^)    &(r  -  r^)   .   (3-7) 
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The  boimdary  conditions  on  <t)  (r)  will  be  the  same  as  those  used  for  the 
thermal  neutron  flux  previously;  the  boundary  conditions  on  q(r^u)  will 
be  identical^  for  u  >  0,  except  that  at  z  =  0^  q(£;'u)  will  be  taken  as 
zero  on  all  boundaries.   This  corresponds  to  assuming  that  <\>{rj-a)   =  0  at 
all  boundaries  for  u  >  0,  which  is  consistent  with  the  thermal  neutron 
flux  assumptions. 

Defining 


lOJ 


2  (u)  +  -7-v 


(3-8) 


Eq.  (3-6)  may  be  written  as 


D(u) 


rffir  ^  ^  ^("'"'  ^  k 


q(r,u) 


2  Jo)    ^ 

=  rw  ^  I  '^"^  ^k  ^^^k^  ^^^  -  ^k)   •  (5-9) 


Multiplying  through  Eq.  (3-9)  by  the  integrating  factor 


u 


exp 


/  p(u',co)  du' 


,  Eq.  (3-9)  may  be  written  as 
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/  p(u'^cj)du' 


u 


ou 


/  P(u',w)du' 


q(r,u)  e 


u 


/  P(u',w)du' 


=  e 


z^(o) 


5(u)  7^   "'s^-k'*  ^^- 


^k) 


(5.10) 


k 


In  order  to  simplify  the  notation_,  define  the  Fermi  age,  i,    from 


dr  _  D(u) 
du   i    2^(u) 


with  t(u  =  0)  =  0 


(5.11) 


and  also  define  a  freq^uency -dependent  resonance  escape  probability, 
p(u,w),  as 


u 


P(u',cj)du' 


p(u,cj)  =  e  ° 


(5-12) 


Notice  that  the  above  definition  is  similar  to  the  usiial  resonance  escape 
probability  except  that  an  additional  l/v  absorber  has  been  added  by  the 
oscillating  source. 

With  the  two  definitions  [Eq.s.  (5-ll)  and  (5-12)]  Eq.  (5-10)  "be- 
comes 


—7 r  ^  q(r.T)  -  5— 


q(£^T) 
_p(t,w) 


T     2.(0) 
J.      .3 

p(t,w)  TToJ  "i 


5(x)  7,  *  (r  )  5(r  - 


k  s  -k 


i.) 


(3-13) 
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The  partial  with  respect  to  the  Fermi  age  will  be  eliminated  by- 
using  the  Laplace  transform.   Defining 


CO 

'(l,s)s      e     r  ,A  dT  ,  (3.11^) 


P(f,w) 


Eq.  (5.15)  becomes 


q(r,Oj     ^^\^)         \' 
^   0(r,s)  -  s  0(r,s)  4-  -^^-^^  .  _       ^  ^  ,^  ,Jr^)  5(r  -  r^)   . 

k  (3.15) 


The  term  q(r^O)  is  zero  since  q(r,0)  was  incorporated  into  the  original 
equation  and  is^  of  course,  the  negative  of  the  right-hand  side  of  Eq. 
(3-15) •   If  it  had  been  chosen  that  this  soirrce  be  introduced  as  a  boundary 
condition,  it  could  have  been  omitted  from  the  original  equation  and  have 
been  introduced  naturally  at  this  point. 

The  obvious  step  at  this  point,  consistent  with  previous  methods, 
is  to  expand  6   in  the  y  eigenfunctions  given  in  Eq.  (2.2),  which  were 
determined  by  Eq.  (2.1)  and  the  boundary  conditions  of  symmetry  and 
•»]/,  (+b)  =  0.   The  required  expansions  are 

0(£,s)  =  ^  0_g(x,z,s)  ^l;_g(y)  (3.I6) 


and 


*s^-^  =  7   %(^^^)  ^^^y)    •  (5-17) 
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For  future  use,  it  should  be  noted  that  Eq.  (3.I6)  implies  that 


q(£,t) 


=  )  Q^(x,z,t)  t^gCy); 


(5.18) 


Eq.  (5-l6)  may,  in  fact,  be  obtained  from  Eq.  (5.I8)  by  Laplace -transform- 
ing both  sides  of  Eq.  (3.I8),  where  Q  is  the  inverse  transform  of  0.. 

Inserting  Eqs.  (5.I6)  and  (j.l?)  into  Eq.  (5.I5)  and  operating  with 


D 

/  t.(y)  dy  gives 


-b 


^.^-(=.^) 


ax^     az2 


i' 


0_g(x,z,s) 


zjo) 


(5.19) 


where 


B^  = 


(2i  -  l)rt 
2b 


(5.20) 


Now  0^(x,z,s)  is  expanded  in  the  x  eigenfunctions  of  Eq.  (2.l4), 
defined  in  Eq.  (2.15) : 


e^(x,z,s)  = 


y  {0,(z,s)  ^cos-^2j-l)nx^g,(    )  l_3i^jjix| 


(3.21) 


Inserting  Eq.  (3-21)  into  Eq.  (3.I9)  results  in 
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v^/ 


(2.1    -   l)n 
2a 


-,2 


-^^S^^I 


1-  "i  i  ii 

a 


sm 


IJTX 


1    1_  COS   ^^J  :  ^^"^  +     )     i     ^  <.-,•.  i22£ 


3z^       Va 


2a 


sm 


az=      '-        =" 


TTX 


Va- 


^s^°)  . 

2^T)    ^    7^6(x   -x^)    5(z    -   .^)    *_g(x,z) 


(5.22) 


a 


Operating  on  Eq.    (5-22)    with      /      —  cos  -^ — dx  gives 


-  e,(z,s) 


(2j    -   l)fl 
2a 


2        c320.(z,s) 

+  — ^ (s  +  B^)    e.(z,s) 


az^ 


^i'      j 


23(0)  ^  (2j    -  Drtx^ 

zToy  ^  L  'k  ^^^  -  \)  7  ^°^  — 2i —  '/V^^ 

t  ^  Va 


(3.23) 


a 
Similarly,  operating  on  Eq.  (3.22)  with   /  —  sin  dx  gives 


e:(z,s) 


2    ^2 


a 


^ (s  +  b2)  0'(z,s) 


2g(0)      ^  ^         JTTX^ 

k 


a 


(3.2U) 
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Defining 


^i^  s  s  +  B?  + 


(2j  -  1)^' 
2a 


(5.25) 


and 


^2  =  s  +  Bf  + 


-,2 


a 


(5.26) 


Eqs.  (5.25)  and  {^.2h)   may  be  rewritten  as 


—J ^^  e.(z,s) 


2^(0)  ^  ^  (2j    -   l)«x 

ZToy  ^     i    \  ^(^    -   \)    -  =°^  2i "oW:^)         (5-27) 


\g^"k^ 


and 


— ^ I  e:(z,s) 


2J0) 

1705- 


I  \  "(^  -  \)  7=  -"  ^  %(v^'    •  (5-==8) 


va 


The  variable  s  has  now  been  hidden  in  n  and  i;    so  the  above  two  equa- 
tions are  effectively  one  dimensional.   It  may  also  be  noted  that  Eqs.  (5.27) 
and  (5.28)  are  Green's  fixnction  eqiaations.   Consider  the  boundary  conditions 
which  are  appropriate  for  these  equations.  Aside  from  the  usual  conditions 
of  continuity  at  z  =  z  and  the  discontinuity  in  the  first  derivative  at 
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that  point,  the  appropriate  conditions  must  be  derived  from  the  original 
z  -axis  condition  on  q(r, t),  which  are 

^(1,-^)1        =  ^i^  ^(Z.-^)  =  0^  Wo  .  (5.29) 

z=o 

Using  these  relations  along  with  the  expansions  that  have  been  utilized, 
Eqs.  (3.16),  (3.18),  and  (3'2l),  results  in  the  following  boundary  condi- 
tions for  Eqs.  (3.27)  and  (3.28): 

e.(z,s)|  :=   ^^™  0.(z,s)  =  0  (3.30) 

J   •*    '        Z-*0O   J        ' 

^  z=o         ^ 

and 

e:(z,s)|   =  ^^™  e'.(z,s)  =  o  .  (3.31) 

"^      z=o         "^ 

Except  for  the  fact  that  the  function  itself,  rather  than  the  first 
derivative,  is  zero  at  z  =  0,  we  now  have  a  problem  identical  to  the  one- 
dimensional  Green's  functions  of  Chapters  I  and  II.  By  following  procedures 
identical  to  those  used  in  Chapter  I,  it  is  easily  shown  that  the  solution 
to 

£!^ifif2l.K^  e(z,zo)  =  -5(z  -  zo)  (3.52) 

dz^ 

with  boundary  conditions 


e(z,Zo)|    =  ^^  0(z,zo)  =  0  (5.55) 

z=o 


IS 


i+0 


-Kz 
ni          \        s    sinh  K  Zn       ^ 
e(z,Zo)  =  j^ ^  :,    z  >  Zo   , 


-Kzo 

e sinh  K  z 

K 


(3.5^^) 


,    Z  <  Zq 


Comparing  Eqs.  (2.25)  and  (2.28)  with  Eqs.  (5.52)  and  (5.5U)  shows  that  the 
change  in  boundary  condition  has  simply  resiolted  in  the  cosh  being  replaced 
by  the  sinh.  By  writing  the  sinh  as  the  difference  of  two  exponentials  in 
Eq-  (3-5^)^  it  is  easily  seen  that  the  two  solutions  may  be  combined  into 
a  form  which  will  prove  to  be  useful  later: 


-K  z-Zr 


e(z,zo) 


2K 


-K(z+Zo) 


2K 


(3.35) 


Using  this  result,  the  solutions  to  Eqs.  (3.27)  and  (3.28)  may  now 
be  immediately  written  down,  being 


2  (0) 

s^  '      ri_ 


^j<^'^'  =  ^  \ 


-\i.\7,--L^  -^i(z+Z  ) 

e         -  e 


7,  —  cos 


(2j  -  l)nx^ 


2a 


%(VV      (3.36) 


and 


e'(z,s) 

J 


s     H— 
270]"  2i 


•^|z-Zvl 


-|(z+Zj^) 


1    .  j"^ 

v/a 


%(v\) 


(3.57) 


i^l 


Inserting  Eqs.  (5-36)  and  (5-57)  into  Eq.  (5 -21)  gives 


'k 


-(J,  I  z-z 


k' 


-M.(z+z^) 


J  k 


(2.1  -  i)nx       l!Li^l!!k, ,   , 


-I  z-z,      -|(z+Zi  ) 

e        -e         .jjTX.k,/     Ni 

+  sm  ^ —  sm  <p„{x,  ,z,  J  r 

a       a    i   k-"  k  I 


(5-58) 


The  next  step  is  to  inverse -Laplace -transform  Eq.  (5o8)^  obtaining 
the  expansion  coefficient  in  Eq.  (5.18)^  which  is 


(5.39) 


Rewriting  the  definitions  of  \i   and  |  as 


H^  =  s  +  B^^   and  ^^  =    s  +   B]^   , 


(3.^0) 


where 


B' 


■ij 


(2i  -  l)n 
2h 


(2j  -  l):t 
2a 


and  B'^ 


(2i  -  1)]T 

2b 


-,2 


a 


(5.i+i) 


it  is  evident  that  each  teim  to  be  inverse-transformed  is  of  the  form 


ys~T~a 


(3.i+2) 


k2 


which  has  the  inverse 


-Qt 


e -{t."^  jhi) 


for  z-^  >  0 


il>M) 


Using  Eq.  (5-^l),  the  inverse  transform  of  Eq.  (5- 38)^  Q«(^^2^''^)^  is 


2^(0) 


^t^^^  2a  y^ 


e  -  e 


J  k 


e      cos 


:^j"     (2.1  -  D^x     ^^J  -  ^^"^  ,   -^i!"^   .   inx   •   "^""^ 

^  nna  -^  V  — i —  cos +  e      sm  ^ —  sm  

2a  2a  a       a 


^^V\) 


(5.H) 


Returning  now  to  the  thermal -neutron  diffusion  equation,  insert  Eqs. 
(5.17)  and  (5.18)  [with  t  replacing  t  in  Eq.  (5.I8)]  into  Eq.  (5.7)  and 


operate  with   /  \!r.(y)  dy  to  obtain 


-b 


.  D  ('t  +  ^  V  i^D  B^  +  Z   +  i^ 


"  'ax^   Sz^  ^    ^  "  ^    ^"   ^ 


*,(x,z) 


s  ^  -" 


=  p(t^5^u))  Q^(x,z,t^)  -  y   7j^  <t>_g(x,z)  5(x  -  x^^)  6(z  -  z^)   .   (3-^5) 


Using  the  expansions 


.,(x,.)=^^[,.,(z)cos(^i^^...^(.) 


ya 


sm 


ijtx 


(5.i^6) 


1^5 


and 


va   . 


i.(z,T^)  COS  ^^'^  ~^^^^^  +  Qj^^^^s^   "^^  ^ 


ItX 


(3.^7) 


in  Eq.  (5-^5)  and  operating  with   /  —  cos   '^   Z — - —  dx  gives 


va. 


2a 


D 


(2,1  -  1)« 
2a 


*.,(z)  -  D 


d2<t)..(z) 
Ji 

s    -,2 


D  B^  +  2   +  — 
s     £  as   V 


^.,(z) 


=  Vi^^,'^)    Q.(-.-3)-^  ^  7^  COS 


(2j  -  l)jtx. 


k 


■i'v^)  ^(^  -  \)  -('-w) 


k 


a 


Similarly  expanding  and  operating  with   /   —  sin  ^ —  dx  gives 


-a 


D 


11 
a 


d^o'  (z) 

<^''  (z)  -  D  —^ + 

J^       "    dz2 


D  B^  +  2   +  — 
s     £  as   V 


*m(^' 


p{t^,u)  Q:(z,-t  ) 


va 


J«Xn, 


7.  sm  ■  fl\"i 

k      a    i   k 


(x  ,z)  5(z 


\) 


(j-w 


Defining 


K^  =  B^  .  + 


iio 
2   +  — 
as   V 
s 


^j     D^ 


(5-50) 


and 


kk 


2    _ 


K'^  = 


Z       +  — 
as       V 


(5.51) 


Eq.s.    (5.^8)    and   (5-^9)   may  be  written  as 


d^O.-.Cz) 


dz 


-P(t    ,w) 


^-■^V^'=        D. 


*J<^'%' 


Dv/a 


(2j    -   l)itx^ 


^k  ^°"  2a 


%K^")    &(z    -  z^) 


(3.52) 


and 


d^*!„(z) 

ai K'2  *!,(z)  = 

dz2  J^ 


pC^,^'^) 


D  ^j(^'S' 


(5.55) 


Comparison  of  the   expansion^    Eq.    (5-^7),    with  Eq.    (5-^^)    shows   that 


Qj^^'-^s^    =   2J0T 


Z   (0)  "^j% 

s  Tie        '^ 


(z-z^)= 


(z  +  2 


k' 


t^^^        2/ 


itax 


-  e 


s         k 


(2j    -  l)nx. 


cos 


k 


2a  %^V\^ 


(5-5^^) 


and 


J+5 


Q!(z,T^) 
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(5.55) 


Since  Eq.  (5-52)  along  with  Eq.  (5*5^)  is  so  similar  to  Eqs.  (5.55)  and 
(5.55)^  attention  will  be  restricted  to  the  former  set;  then  the  latter 
set  will  he  handled  by  comparison. 

Equations  (5-52)  and  (5.5^)  may  be  combined  as 
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(5.58) 
(5.59) 

Equation  (5-56)  is  virtually  identical  to  Eq.  (I.50)  and  has  es- 
sentially the  same  boundary  conditions,  so  that  the  Green's  function  of 
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Chapter  I  may  be  used  to  convert  Eq.  (3- 56)  to  an  integral  equation.  The 
appropriate  Green's  function  is  determined  by 


S^G.Cz,!) 

J K^  G,(z,^)  =  -  6(z  -  I) 


(5.60) 


and  the  usiial  homogeneous  boundary  conditions.  The  solution  to  Eq.  (5-6o), 
from  previous  results,  is 


G.(z,0 
J 
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(5.61) 


As  before,  the  integral  equation  for  <l>..(z)  is  foimd  by  forming 

00  CO 

/     d|  Gj(|,z)    [Eq.    (5-56)]    -    J      ^1   ^j^(0    [Eq.    (5-60)],    after  reversing 
o  o 

z  and  I   in  Eqs.    (5.56)    and  (3.6O).      This  resvilt  is 
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An  integration  by  parts  may  be  performed  on  the  integrals  containing  the 
derivatives  and  the  integral  containing  H   may  be  completed^  resulting  in 
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(3.65) 


Inserting  the  boundary  conditions  into  the  last  terra  gives 
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Inserting  A   and  H.,  and  assuming  that  z^,  <  z  <  z^,_^^^  Eq.  (5-6^)  becomes 


M' 


M'+l' 


».,(z) 


el^  "♦,W'^) 


M' 


K      d| 


|=o 


KD^a 


-Kz 


cosh  Kz,  7,  cos 
k  k 


(2j  -  l)nx^ 


2a 


k=l 


M 


<l5„(x,  ,z,  )  +  cosh  Kz 
^   k'  k 


-Kz         (2j  -  l)rtx. 


k 


^k  ^°"  2^ 


k 


:^V"k^ 


k-M'+l 


ZjO)  n  p(t  ,03)  e  ^J  " 


Z^(0)  2  K  Dg/^ 


-Kz 


d|  cosh  K|  2^  Fj^d.Zj^) 


CO 


-K^ 


+  cosh  Kz   /   d^  e 

z  k 


I  ^Jlc(5'\' 


(3-65) 
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It  seems  appropriate  at  this  point,  before  continuing  into  the 
mathematical  complexities,  to  simply  the  notation  in  Eq.  (3.65)  and  try- 
to  -understand  its  physical  significance.   To  this  end  call  the  first  term 
in  Eq.  (5 -65)  H..,  since  it  represents  the  contribution  -which  would  be 
present  in  the  homogeneous  case.  Also,  redefine  Eq.  (3.58)  aJi<i  (3-59)  a-s 


F   (z,z^)  =y^Mj^(z)  W   %(vzj^) 


(3.66) 


and 


«J.^^) 


D  V  %K'"k) 


(5.67) 
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W.T  =  —^   cos 
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■/a 


(5.69) 


W   is  simply  a  -weighting  factor  depending  on  the  distance  from  the  center 
Jk 

line  of  the  foil  location.  M,  (z)  is  the  slowing-down  kernel  for  the  parti- 
cular geometry  chosen,  having  an  obvious  resemblance  to  the  infinite -medium 
Fermi -age  kernel.  Noticing  the  factor  ij  Jl~ ,    it  can  be  seen  that  the 
slowing-down  kernel  obtained  for  the  present  geometry  resembles  the  plane 
soirrce  infinite -medixim  slowing-down  kernel  most  closely.   This  is  to  be 
expected,  since  the  problem  was  reduced  to  a  semi -infinite,  one -dimensional 
case. 


1^9 

Inserting  Eqs.  (3'66)  and  (5-67)  into  Eq.  (5.65)  and  using  the 
symmetric  part  of  Eq.  (3- ^6)  and  Eq-  (5 -61)  results  in 


^,(x,z)  =  ^  W.(x)  H.,  -  ^)^  _il_w.^G.(z^,z)  */x^,z^) 


00 

I  I^^'ji  ^J^"^  ^Jk  ^^v^)  /  ^^  V^^  ^j^V^)    ^  ^5.70) 


where 

W.(x)  =  l,cos  (2i^^   .  (5.71) 


It  should  now  be  noted  that  the  ith  mode  of  the  flux  distribution  consists 
of  the  following  parts:   (a)  the  homogeneous  term  /   H^^  which  represents 

the  propagation  of  a  neutron  wave  in  a  homogeneous  system^  (b)  the  absorb- 
ing part^  the  second  term^  which  accounts  for  the  heterogeneous  absorption 
of  the  plates,  and  (c)  the  last  term,  which  accounts  for  the  production  of 
fast  neutrons  by  fissions  in  the  plates,  and  their  subsequent  slowing  down 
into  the  thermal  group. 

Returning  to  the  mathematical  manipulations,  it  now  remains  to 
evaluate  the  two  integrals  in  Eq.  (5'65),  which  are  complicated  by  the 
fact  that  K  is  a  complex  quantity.   The  first  integral  may  be  written  as 


z 

r 

\~^ 

r "  ' 

Ut 

/    d^    cosh  K^ 

) 

t 

s 

0 

k 
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where 


(3.75) 


Expressing  the  cosh  as  the  sum  of  exponentials  and  performing  the  squares 
in  the  exponentials,  Eq.  (5-72)  may  be  reformulated  as 
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Multiplying  and  combining  the  various  exponentials  leads  to 
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By  completing  the  squares  in  each  of  the  exponentials,  Eq.  (5-75)  may  be 
rewritten  as 
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(3.76) 


Each  of  the  four  integrals  in  Eq.  (5- 76)  is  of  the  form 
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(5.77) 


which  becomes,  upon  changing  the  variable  of  integration  to 


w  =  -^I:+^AT  a  , 
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(3.78) 
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Recognizing  the  two  integrals  as  error  functions^    Eq.    (5'79)   becomes 
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where  the  error  fimction  is  defined  as 
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(5-81) 


The  integrals  of  Eq.  ('J-TS)  now  become 
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The   second  integral  in  Eq.    (5.65)    is 
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which  can  be  handled  in  a  similar  manner  to  give 


(3.82) 


(3-85) 
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If  the  squares  in  the  factors  exp 
{^.8k)    are  performed,  a  factor  exp 
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in  Eqs.  (3.82)  and 
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will  res-Qlt  which  will  cancel 
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The  two  remaining  factors  are  exp 
,    which  can  be  taken  out  of  the   summation  and  combined  with 
the  factor  exp     -  B^ .    t 


Now  the  expressions  for  the  two  integrals,    Eqs.    (5.82)    and   (5-8^)^ 
may  be  inserted  into  Eq.    (5.65)    to  give  the  complete  expression  for  <t)      (z) 
Thus: 
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(continued  on  following  page) 
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(5.85) 


Equation  (5.85)  with  Eq.  (5. ^6)  used  for  <t>  (x,  ,z  )  is  the  complete 
solution  to  Eq.  (5.52).   <t>I  (z)^  the  solution  to  Eq.  (5-53)^  is  similar  to 
Eq.  (5.85)  except  that  K'  replaces  K,    B' .  replaces  B.  .,  and  sin(jitx,/a) 


replaces  cos 


(2j  -  l)rtx^ 
2a 


wherever  they  appear.   Equation  (3.85)  and  the 


corresponding  equation  for  Ol  .(z)  are  inserted  into  Eq.  {J).h6),    which  in 
turn  is  inserted  into  Eq.  (j-l?)  to  give  the  complete  expression  for  <t>(r)^ 
the  thermal  neutron  flux  at  any  point  in  the  medium. 

It  is  instructive  to  see  what  happens  to  Eq.  (5-85)  when  the  age- 
dif fusion  theory  is  reduced  to  one -group  diffusion.   If  no  errors  have  been 
made^  Eq.  (5.85)  should  reduce  to  essentially  the  same  form  as  Eq.  (I.66). 
To  do  this,  let  t  ->  0  and  set  p(t  ,oj)  =  1.   As  t  -*  0,  the  arguments  of 


all  the  error  functions  will  approach  the  form  or  +  


Oh- 


2^/t 


2Jx 


viously  all  the  arguments  will  approach  +00  on  the  real  axis,  depending  in 

some  cases  on  whether  z  <  z,  or  z  >  z,  .   The  error  function  for  +00  is 

k        k  — 

simply  +1,  and  so  the  last  term  in  Eq.  <3.85)  becomes 
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(5.86) 


With  a  little  further  manipiilation  Eq.  (3.86)  can  be  further  reduced  to 
the  form 
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k=M'+l 


Insertion  of  Eq.  (5.87)  into  Eq.  (3'85)  reduces  the  latter  equation  to 


57 


-Kz   d*   (I) 
<l>..(z)  =  ^ 


j^' 


K 


d| 


l=o 


23(0) 

z^oy  ^  - 1 


1     e 


D^V^ 


-Kz 


K 


M' 


(2j  -  l)nx. 


cosh  Kz,  7,  cos — 

k  k        2a 


*/v\^  ^ 


cosh  Kz 
K 


k=l 


M 


-Kz 


k 


7,  cos 
k 


(2j  -  l)jtXj^ 


2a 


%^V\^ 


k^M'+l 
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The  resemblance  of  Eq.  (3.88)  to  Eq.  (I.66)  is  obvious.   Except  for 
the  differences  in  the  shapes  of  the  perturbing  foils  assumed  in  the  two 
derivations^  by  following  a  procedure  similar  to  that  described  at  the 
end  of  Chapter  .II  the  equivalence  of  the  two  expressions  can  easily  be 
shown. 


CHAPTER  rV 


CALCUIATIONAL  RESULTS 


Two  computer  codes  were  written  to  obtain  numerical  results  from 
the  analytical  results  of  the  previous  chapters.   The  simpler  of  the  two 
uses  the  results  of  Chapter  I  with  no  additional  analytical  simplifications, 
The  second  code  uses  the  age-diffusion  model  developed  in  Chapter  III,  but 
several  restrictions  had  to  be  made  because  of  computer  storage  and  numeri- 
cal problems. 

The  simpler  of  the  two  codes  performs  the  procediires  given  from 
Eq.  (1.67)  to  the  end  of  Chapter  I.   Computer  storage  limitations  resulted 
in  several  lijnitations  on  variable  sizes.  M,  the  number  of  foils,  was 
limited  to  6;  using  the  terminology  of  Eq.  (1.6?) ,  L  was  limited  to  10,  Q 
was  limited  to  5,  and  the  product  LM  was  limited  to  50.  The  output  of  the 
code  consisted  of  0^^(z)  for  up  to  I+9  uniformly  spaced  values  of  z  and  the 
summation  of  all  the  flux  modes  for  the  same  z  values.   These  results  were 
given  by  the  code  for  any  number  of  \aniformly  spaced  frequencies  of  the 
driving  source. 

Since  the  one -group,  one -dimensional  computer  code  almost  completely 
filled  the  computer  core,  it  was  obvious  that  the  two-dimensional  age- 
diffusion  res\ilts  could  not  be  coded  in  their  entirety.   Two  major  simpli- 
fications were  made.   One  was  to  restrict  the  foil  positions  and  flux 
calculational  points  to  the  z  axis  (x^  =  x  =  O),  thereby  transforming  the 
problem  to  a  one -dimensional  problem.  This  results  in  the  x  expansion, 
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Eq.  {^.k6),    reducing  to 
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Epilation  (3'T9)  simplifies  somewhat  since  the  factor  cos 

"becomes  unity,  and  since  <t)'.  .(z)  is  not  used  in  Eq.  (^.l),  the  equation  cor- 
responding  to  Eq.  (5-79)  for  '^'.gi'z-)    is  not  needed.   Substitution  of  Eq. 
(^.l),  with  z  replacing  z,    into  Eq.  (3'85)  along  with  replacing  the  cosine 
factors  with  unity  gives  the  basic  equation  the  code  must  handle  as  follows; 
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The  only  other  simplification  that  was  made  \j3.s   only  the  inclusion 
of  the  fundamental  spatial  mode  in  the  source.   Formally^  this  implies  that 
Eqs.  (1.13),  (1.22),  and  (l.2U)  apply  and  that  S   s  0  for  i  >  1  and  m  >  1. 
In  practice,  the  higher  i  values  (x  modes)  could  have  been  easily  included, 
but  to  have  included  the  higher  m  values  (y  modes)  woiold  have  overloaded 
the  computer  core.  This  is  not  an  insurmountable  obstacle,  since  magnetic 
tape  storage  could  have  been  used  but  this  is  a  time-consuming  and  there- 
fore expensive  approach.  Physically,  including  only  the  fundamental  spatial 
mode  seems  justified  since  recent  experimental  work  at  the  University  of 
Florida  indicates  that  a  highly  thermalized  fundamental  mode  is  experi- 
mentally feasible  (lO) . 

Using  Eq.  (U.2)  rather  than  Eq.  (I.66),  the  procedure  outlined  in 
Chapter  I  from  Eq.  (I.66)  to  the  end  of  the  chapter  may  be  used  to  obtain 
the  complete  solution  to  the  problem.   This  will  be  done  later  in  this 
chapter,  but,  first,  attention  will  be  directed  to  two  factors  in  Eq.  (U.2) 
which  have  physical  significance  that  has  not  been  pointed  out  previously. 

If  Eq.  (3.8)  is  inserted  into  Eq.  (3. 12),  with  the  age  to  thermal. 


T  ,  replacing  u,  the  result  is 
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p  s  Z  (u'  )du'        p  s     ^  , 
/    a  .    /       du' 

"  J    |2^(u')    -"'^  J    |2^(u')v(u') 
p(t^,w)  =  e  °  e    °  ,  (h.^) 

Noting  that  the  first  integral  in  Eq.  (^.j)  is  the  usual  resonance  escape 

probahility  to  thermal^  p(t  )  or  p(u  )^  Eq.  (^.j)  may  "be  written  as 

s        s 

-iwT 
p(t  ,0))  =  p(t  )  e    ^   ,  (^.U) 


where 
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or,  in  terms  of  the  Fermi  age, 

T 

s 
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Ohviously,  T  has  the  units  of  time,  and  represents  the  time  required  for 
the  neutron  to  slow  down  to  thermal  energy.   The  factor  exp  [-itoT  ]  then 
represents  the  phase  shift  introduced  by  the  slowing-down  time. 

The  second  factor  in  Eq.  {k.2)    which  can  be  simplified  is 

exp[(K^  -  B^.)t  ].   With  the  use  of  Eq.  (3.5O),  this  immediately  becomes 
zj      s 

icj 
2  +  ^ 

as  v^ 

(K^-Bf.  )t        D     ""s  (k.7) 

£j      s  s 

e      ^    =  e 

2    T  T 

as   s     .     s 


D         V  D 
e   ^    e     ^   ^   .  (U.8) 
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The  term  in  the  first  exponential  is  the  ratio  of  the  age  to  thermal  to 

the  square  of  the  diffusion  length,  t  /l^.   Remembering  that  for  a  point 

s 

source  L  is  one-sixth  the  mean  square  distance  that  a  thermal  neutron 
travels  before  being  absorbed  and  that  t(u)  is  one-sixth  the  mean  square 
distance  that  a  neutron  travels  in  attaining  lethargy  u,  the  ratio  can  be 
seen  to  be  a  measure  of  the  distance  traveled  while  slowing  down  relative 
to  the  distance  traveled  at  thermal  energy.   Typical  values  for  this  ratio 
are  O.15  for  graphite  and  O.OO8  for  E^O. 

Now  returning  to  the  actual  calculational  procedure,  as  before,  the 
expansion,  Eq.  (l+.l),must  be  truncated;  i.e., 
J 

^^  J=l 

Because  of  the  restriction  of  the  source  to  the  fimdamental  mode  only, 
Eq.  (4.9)  simplifies  further,  since  £   can  only  assume  the  value  1. 
Restricting  j  to  1  in  S   does  not,  of  course,  eliminate  the  higher  j  com- 
ponents in  the  flux  since  the  perturbing  foils  generate  higher  x  modes  but 
they  cannot  generate  higher  y  modes  since  the  foils  are  assumed  to  extend 

the  full  length  of  the  medium  in  the  y  direction.   Equation  (5.85)  with 
J 

—   2_,  *jl^^k^  inserted  for  %(\^Zj^)  and  cos  — ^  replaced  by  unity 

Va  j=l 

may  now  be  written  JM  times  for  each  of  the  J  modes  and  for  each  of  the  M 

foils.   This  set  of  equations  may  then  be  solved  simultaneously  for  the  JM 

values  of  '^■ji^)-      Writing  these  values  as  a  vector, 

""^^  ^^l^^l^'  *1i(^2^^---^i(^m)^  W\^'---'j1^\^^      >  (^-10) 
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the  set  of  JU   eq_uations  may  be  put  into  matrix  form,  as  before.   That  is, 


[A]  $  +  [I]  $  =  B  , 


(^•11) 


where  B  is  the  vector  of  source  terms,  [I]  is  the  identity  matrix,  and  [A] 
is  the  matrix  containing  the  Green's  functions  connecting  the  foils.   B 
may  be  written  as 


B^  = 


a  ,  e      ,  . . .  a   e      ,  0  ...  0 


(^.12) 


where 


a 


11 


11   2DK. 


11 


(^.15) 


The  [a]  matrix  is  defined,  as  before,  by  partitioning  it  into  J^  submatrices, 
each  M  by  M.   The  n,p  element  of  the  I,K  submatrix  is  given  by 

— K   z 
^^I,K)n,p  =  ^Il^^p)  ^°^^  ^Il^P  ^   '  "  -  ^1^  Vn,p^   for  n  >  p  ,  (k.lk) 


=  6^^  (z  )  cosh  K^-,z   e 
II  p        II  n 


-h 


1  P 


Wn,p^   forn<p   ,   (^.15) 


where 


i^J   - 


II'  p'    K  Da  ' 


{h.l6) 


2  + 


ICJ 


^1- 


D 


2  D  a  K  , 
s    II 


(^.17) 


6h 


-K_^(z    -z    )    7 
(Vn,p  -  e    ^1     ^     P     ^     ^(Al)    -   (A2)    4-  2(A5) 


-KJz+z)    7 


+  e 


'^     "     ^     2^     {(Al.)    -    (A5)    -  2(A6)| 


+  cosh  K^^z^  7        je   ^-^  P    [(A2)    -   l]    -  e      ^^  ^    [(a1|)    -   U        ,      (^.l8) 


(Al)    =  erf 


z      -  z 

^-^^K,1 

2yr       ^  II 


(^.19) 


(A2)    =   erf 


z      +   z 


(i|.20) 


(A5)    =   erf 


2^7" 
^    s 


+  Vx     K. 


s      II 


(i+.2l) 


(Ai|)   =   erf 


z      -   z 


2yr       ^  ^1 


(U.22) 


(A5)    =   erf 


z      +   z 


2/r       ^ ^i 


(U.25) 


(A6)   =   erf 


2jT 
s 


^^1 


(i+.24) 


The  vector  $  is  obtained  from 


[A]  +    [I] 


B 


(^.25) 
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The  elements  of  $  are  inserted  into  Eq.  (3'85)  with  the  previously  mentioned 
simplifications,  which  is  in  turn  inserted  into  Eq.  (^.9)  "to  give  the  com- 
plete answer  to  <^^{0,z)  .      The  entire  thermal  neutron  flux  is  obtained  by- 
using  Eq.  (3.17),  which  gives 


*  (0,0,z)  =  —  *(0,z)=-^   )   0   (z)   .  (I4.26) 


Since  the  representation  of  the  flux  is  by  expansion  in  eigenfunc- 
tionSj,  an  Important  characteristic  of  the  numerical  calculations  is  the 
rate  of  convergence  with  increasing  numbers  of  the  eigenfunctions,  or 
spatial  modes.   Figures  2  and  5  show  flux  amplitude  versus  distance  along 
the  z  axis  with  1,  2,  k,    and  10  spatial  modes  and  with  the  source  modulated 
at  100  and  1000  cps,  respectively.   The  model  is  the  one-group  representa- 
tion of  Chapter  I  with  the  geometry  shown  in  Fig.  1.   One  cadmium  strip, 
0.0508  X  1  x  71-12  cm,  is  positioned  8.89  cm  from  the  end  of  a  graphite 
assembly  !+il-.^5  x  71-12  cm,  the  dimensions  of  one  of  the  experimental 
facilities  available  at  the  University  of  Florida.   In  addition  to  the 
curves  shown  in  Figs.  2  and  5^  calculations  were  performed  with  6  and  8 
spatial  modes.   At  the  position  of  slowest  convergence,  which  is  in  the 
neighborhood  of  the  cadmium  foil,  the  ratios  of  the  6-  and  8-mode  results 
to  the  10-mode  results  are  approximately  I5  and  6^,  respectively,  for  both 
100-  and  1000-cps  source  frequencies.   Figure  ^4-  shows  the  phase  lag  versus 
z  for  100  cps,  corresponding  to  Fig.  2.   In  this  case  the  change  from  6  to 
10  modes  and  from  8  to  10  modes  is  approximately  5  and  2/0,  respectively. 
Figure  5  is  the  phase  lag  versus  z  for  1000  cps,  corresponding  to  Fig.  3« 
The  convergence  is  faster  for  this  case,  approximately  2.5  and  ifo   for  6  to 
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Fig.  2.   One -Group  Flux  Amplitude  vs  Distance  Along  the  Z  Axis  for 
One  Cadmi-ura  Foil  in  Graphite,  at  100  cps,  with  1,  2,  \,    and  10  Spatial 
Modes. 
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Fig.  3'  One-Group  Flux  Amplitude  vs_  Distance  Along  the  Z  Axis  for 
One  Cadmium  Foil  in  Graphite,  at  1000  cps,  with  1,  2,  h,  and  10  Spatial 
Modes. 
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Fig.  h.      One -Group  Flux  Phase  Lag  vs  Distance  Along  the  Z  Axis  for 
One  Cadmium  Foil  in  Graphite,  at  100  cps,  with  1,  2,  \,    and  10  Spatial 
Mode  s . 
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Fig.  5.  One-Group  Fliix  Phase  Lag  vs  Distance  Along  the  X  Axis  for 
One  Cadmium  Foil  in  Graphite^  at  1000  cps,  with  1,  2.,  h,  and  10  Spatial 
Modes. 
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10  and  8  to  10  modes,  respectively.   CalcixLations  were  limited  to  10  spatial 
modes  because  of  machine  storage  limitations.  Altho-ugh  the  convergence  is 
somewhat  slow,  it  is  probably  adequate  for  most  purposes,  especially  since 
the  problem  used  is  a  rather  severe  test. 

Figixres  6-10  show  the  results  of  calculations  in  the  same  geometry 
as  described  above  except  that  two  perturbing  foils  are  present,  at 
z  =  8.89  and  17.78  cm.   These  foils  are  either  cadmium,  as  used  previously, 
or  U^^^  of  the  same  total  absorption  cross  section  as  the  cadmium  strips. 
The  calculations  with  cadmium  are  one -group  and  with  U^^^  are  age -diffusion. 
Figure  6  shows  fliix  amplitude  as  a  f-unction  of  z  for  one  spatial  mode  cal- 
culations using  cadmium  foils  with  source  frequencies  of  200,  500,  and  8OO 
cps  and  using  iP^^  foils  at  200  cps.   The  unperturbed  flux  at  200  cps  is 
also  plotted  on  Fig.  6.   Figure  7  shows  the  same  data  using  six  spatial 
modes.  The  unperturbed  flux  is  identical  in  Figs.  6  and  7  since  the 
assembly  is  driven  with  fundamental  mode  only  in  each  case.   Figure  8 
presents  the  phase  angles  as  a  function  of  z  corresponding  to  the  amplitudes 
shown  in  Figs.  6  and  7-   The  one-group,  six-spatial  mode  calculations  shown 
in  Figs.  7  and  8  are  replotted  as  a  ftmction  of  soiirce  driving  frequency  in 
Figs.  9  and  10.   Figure  9  presents  flux  amplitude  versus  frequency  for  z 
values  of  2,  8,  8.89,  I5,  and  50  cm.   Figure  10  presents  phase  lag  for  the 
same  conditions. 

A  second  facility  that  is  available  at  the  University  of  Florida  for 
measvirements  compatible  with  the  geometry  restrictions  of  these  calculations 
is  a  water,  or  DgO,  tank  which  is  120.97  x  15O.17  cm  in  the  transverse 
directions.   Figure  11  is  a  plot  of  flux  magnitude  versus  z  at  200  cps 
for  an  age-diffusion  calc\alation  in  E^O  with  two  l-in.-diam  natural  uranivrai 
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Fig.  6.   Flux  Amplitude  vs  Distance  Along  the  Z  Axis  for  Two  Cadmium 
or  Uranium  Foils  in  Graphite,  at  200,  500,  and  800  cps,  with  One  Spatial 
Mode. 
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Fig.  7.   Flux  Amplitude  vs  Distance  Along  the  Z  Axis  for  Two  Cadmium 
or  Uranium  Foils  in  Graphite,  at  200,  500,  and  800  cps,  with  Six  Spatial 
Modes. 
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Fig.  8.   Flux  Phase  Lag  vs_  Distance  Along  the  Z  Axis  for  Two  Cadmium 
or  Uranium  Foils  in  Graphite,  at  200,  500,  and  800  cps,  with  One  and  Six 
Spatial  Modes. 
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Fig.  9-   One -Group  Flux  Amplitude  vs  Soirrce  Frequency  for  Two  Cadmium 
Foils  in  Graphite,  at  z  =  2,  8,  8.89,  1'?,    and  50  cm,  with  Six  Spatial 
Mode  s . 


75 


4.0 


3.5 


3.0 


2.5 


o 

o 

< 


2.0 


< 
I 
a. 

X 

3 


1.5 


1.0 


0.5 


./ 

/ 

/z  =  30  cm 

/ 

^/--^ 

cm 

> 

X 

8.89  cm^ 
^^'''^^,-- 8  cm 

/ 

/^ 

=^ 

'^^"^'"^ 

2  cm 

M 

€^ 

200 


400  600 

SOURCE  FREQUENCY  (cps) 


800 


1000 


Fig.  10.   One -Group  Flux  Phase  Lag  vs_  Source  Frequency  for  Two  Cadmium 
Foils  in  Graphite^  at  z  =  2^  8^  8.89^  15^  and  50  cm,  with  Six  Spatial 
Modes. 
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rods  150.17  cm  long  located  at  z  =  10. 998  a.-d  21. 996  cm.   The  obvious 
characteristic  of  the  calciilation  is  that  the  natural  uranium  rods  do  not 
offer  sufficient  perturbation  to  produce  dramatic  effects.  An  alternate 
and  feasible  perturbation  could  be  obtained  by  using  single  plates  from 
fuel  elements  of  the  MTR  type.   One  of  these  plates  has  approximately  the 
same  total  absorption  cross  section  as  the  natural  uranium  rods.  Therefore, 
the  perturbation  would  be  approximately  the  same.   In  order  to  produce  a 
significant  perturbation  with  fuel,  fairly  large  U^^^  rods  would  have  to 
be  used. 

Since  the  calculation  shown  in  Fig.  11  has  small  perturbations,  it 
presents  a  good  opportunity  to  show  the  relative  contributions  of  the 
unperturbed,  thermal -neutron  diffusion,  and  slowing-down  contributions  to 
the  final  answer.   These  are  shown  in  Fig.  12,  where  the  three  curves  cor- 
respond to  the  three  terms  in  Eq.  (5.85).   Since  it  is  a  six-spatial  mode 
calculation,  the  thermal -neutron  diffusion  and  slowing-down  curves  are  sums 
of  the  second  and  third  terms  in  Eq.  (5.85). 
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Fig.  12,   Components  of  Age-Diffusion  Flijx  Amplitude  vs_  Distance 
Along  the  Z  Axis  for  Two  Uranium  Foils  in  E^;0,  at  200  cps^  with  Six 
Spatial  Modes. 


CHAPTER  V 

SUMMARY  Am)  CONCLUSIONS 

The  goal  of  this  dissertation  has  been  to  predict,  as  accurately 
as  possible,  the  behavior  of  neutron  waves  in  heterogeneous  multiplying 
and  nonmultiplying  media.   The  primary  approach  has  been  similar  to  the 
heterogeneous  reactor  theoiy  of  Feinberg  and  Galanin  in  that  the  continu- 
ous slowing-down  model  has  been  ass-umed  to  describe  the  behavior  of  the 
neutrons  above  thermal  energy  and  diffusion  theory  has  been  ass-uoned  to 
describe  the  behavior  of  the  thermal  energy  neutrons.   Rather  than  assume 
that  the  infinite -medium  slowing-down  and  diffusion  kernels  are  applicable, 
a  more  basic  approach  was  chosen  which  resulted  in  slowing-down  and  diffu- 
sion kernels  for  the  particular  finite  geometry  that  was  assumed. 

The  preceding  chapter  showed  some  results  of  calculations  using 
this  theory.  In  this  chapter,  it  will  be  shown  how  this  theory  may  be 
applied  and  extended. 

Application  to  Critical  Reactors 

As  stated  previously,  the  theoretical  developments  of  this  work  are, 
in  several  ways,  an  extension  of  the  Feinberg-Galanin  heterogeneous  reactor 
theory.   Aside  from  the  obvious  feature  of  predicting  the  behavior  of  neu- 
tron waves  in  subcritical  assemblies,  this  work  takes  into  account  the 
finiteness  of  the  assembly.   The  effect  that  this  has  is  best  shown  by 
applying  the  theory  to  a  critical  reactor,  which  can  be  done  by  dropping 
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the  source  terms  and  changing  the  boimdary  conditions  slightly.   The  ab- 
sorptive Green's  function  defined  by  Eq.  (5.6O)  will  have  the  same  boundary 
conditions  as  before  except  that  now  G.(z^|)  itself,  rather  than  the  deriv- 

ative  of  G.(z,|),  will  be  required  to  go  to  zero  at  z  =  0.   This  change 
J 

in  boundary  condition  simply  results  in  G.(z,|)  changing  to 

J 

G.(z,|)  =  £^^-fl^LM     ^        ,>^   ,  (5.1) 

^  e phj^  ^   ^^^   ^  (5^2) 


With  the  new  Green's  function,  Eq.  (5- 63)  applies  to  the  critical  assembly 
with  the  added  simplification  that  all  the  bo^undary  conditions  given  by 
the  last  terra  in  Eq.  (5-65)  vanish.   Equation  (5-65)  now  becomes 

00 

*.,(z)  =  -l^O.i.^,^)    H.^(z)  .A^.  ^   /  d^  G.(^,z)  F.^(|,z^)  ,      (5.5) 
k  k   o 


where  G.(|,z)  is  now  given  by  Eqs.  (5-l)  and  (5-2).   Remembering  that  the 
J 

factors  H.,  and  F.,  both  contain 
jk      jk 


it  can  be  seen  that,  for  k  -   1,2, ...,M,  and  j  =  1,2, ...,J,  JM  homogeneous 

equations  are  available  for  each  value  of  £.      For  the  system  to  be  critical 

the  determinant  of  the  coefficients  of  the  *.  .(z,  )  must  be  zero.   This  of 

j£     k' 

course  gives  an  equation  which  can  be  solved  for  t^  (contained  in  A..)  or 
whatever  other  eigenvalue  is  chosen  to  determine  criticality. 
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The  first  difference  noticed  in  comparison  with  the  Feint erg-Galanin 
theory  seems  to  be  simply  additional  complexity  since  the  Feinberg-Galanin 
theory  only  requires  the  solution  of  an  M  by  M  determinant.  The  improve- 
ments are  the  improved  kernels  contained  in  Eq.  (5-5)  •   As  pointed  out  in 

Chapter  III^  A  .  and  F.,  contain  the  finite -medium  slowing-down  kernel 
^3  Jk 

rather  than  the  infinite -medium  kernel  used  by  Feinberg-Galanin.   The 
Green's  function  G.  also  is  an  improvement,  being  the  finite -medium  diffu- 
sion  kernel  for  the  geometry  used  in  this  work,  rather  than  the  infinite- 
medium  diffusion  kernel  used  by  Feinberg-Galanin. 

Conclusions 

The  theoretical  developments  of  this  dissertation  probably  represent 
a  model  which  is  as  complex  and  physically  realistic  as  is  practical  at  this 
time.  Further  improvements,  in  order  to  be  practical  for  machine  calcula- 
tions, would  certainly  be  desirable  if  they  simplified  the  mathematics  or 
if  they  eliminated  remaining  physically  -undesirable  assumptions,  or  both. 
If  such  improvements  are  not  easily  attainable,  the  primary  question  remain- 
ing is  how  to  use  the  theory  best  to  direct  future  experimental  work  in  this 
area. 

Several  approaches  are  possible.   One  is  simply  to  continue  full 
calculations  of  the  thermal  flux  in  various  configurations,  such  as  those 
presented  in  the  preceding  chapter.   This  is  a  rather  brute  force  technique 
but  at  least  has  the  advantage  that  calculations  are  cheaper  and  much  more 
easily  accomplished  than  experiments  since  the  machine  codes  are  available. 
Some  possibly  more  satisfying  approaches  can  be  s-uggested  by  re -examination 
of  some  of  the  results  of  Chapter  III.   To  this  end,  write  Eq.  (3.6U)  as 
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j,^  (2j    -  l)jtx 


k  ^ 


^    7  (2j    -  l)nv  A^ 


(5.5) 


vhere 

-Kz     d0  . , (I ) 


H.„  = 


e i£ 


j£  K  d| 


(5.6) 
|=o 


-B^.T 


and 


^ij   ^  2 JO)  ^2  D  ^5-7) 


and  G     and  R    are  defined  by  Eqs.    (5.6I)    and   (5.68)^    respectively.      A  cor- 
responding equation  exists  for  <I>I    (z);    thus 

0-^ 


1  V         ^v     j^^v 


k      ^ 


-  e;.  ^  ^  sin  ^  .^(x^,z^)  /  d|  V^)  G'(z^,|)   ,       (5.8) 


k 


Va 


where  G'  H'  ,    and  E'   are  obtained  by  replacing  Bf .  with  Bl^  and  K^  with 
K'^  where  they  appear  in  Eqs.  (5.61),  (5-6),  and  (5.7),  respectively. 
Inserting  Eqs.  (5.5)  and  (5.8)  into  Eq.  (3.^6)  and  in  turn  into  Eq.  (5-17) 
results  in 
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,        .           \    1               (2i    -   l)rty    V    J    1               (2j    -   l)Ttx 
(r,aj)   =      )   —  cos  -^ ;::: — '—^     >      i   —  cos  ^^—^ 


£ 


2b 


H. 


2a  i£ 


1        .      jrtx  „,  1 

+  —  sm  ^-i —  H '.  .    -  - 
a        j^        D 


Va 


.°J<^'^1.'    'iV. 


J      ■*   k       jk 


%^V^k^  -^  ^ij  I  ^jk  ^(v\) 


00  _  oo 


k 


(5.9) 


where 


_  ^k  ^_    (2j    -   l)nx 


C  .,    =  —  cos 

jk       a 


(2j      -     l)lTX^ 


2a  "°"  2a 


(5.10) 


and 


7-, 


S.,    =  —  sm  ^ —  sm 
jk       a  a  a 


j«x^ 


(5.11) 


Experimental  measiirements   of  thermal  neutron  flux  as   a  function  of  position 
and  frequency  represent  the   quantity  calculated  in  Eq.    (5.9)-      These   data 
may  be   operated  on  with 


dx  cos 


-a 


(2j    -   l)nx      r  {2£    -  l)ity 

2a  J    ^^  ^°"  2b 

-b 


(5.12) 


or  with 

a  b 

/    dx   sm  ^ —   /    dy  cos  -^^ _ 

-a  -b 


(5.15) 


Qk 

to  yield  experimental  measurements  of  the  flux  components  calculated, 
respectively,  in  Eqs.  (5.5)  and  (5-8) •   This  entire  procedure  may  be  per- 
formed experimentally  with  the  moderator  alone,  with  fuel  plates  or  rods 
and  with  purely  absorbing  plates  or  rods  of  the  same  macroscopic  absorbing 
cross  section  as  the  fuel  assemblies.   The  technique  of  replacing  fuel  as- 
semblies with  comparable  pure  absorbers  has  been  utilized  previously  to  mea- 
sure effective  delayed  neutron  fractions  (ll) .   The  six  terms  in  Eqs.  (5-5) 
and  (5.8)  may  now  be  separated.   The  homogeneous  experiment  yields  H   and 
H'  .   Subtracting  the  homogeneous  term  from  the  pure-absorber  data  yields 
the  absorption  terms,  the  second  terms  in  Eqs.  (5-5)  and  (5-8).   Similarly, 
the  difference  between  the  pure  absorber  and  the  fuel  assembly  data  yields 
the  last  terms  which  contain  the  slowing-down  kernels. 

Since  the  above  experimental  data  would  presumably  be  obtained  by 
cadmium-difference  measurements  with  a  neutron  detector  which  has  approxi- 
mately a  l/v  response,  it  is  also  of  interest  to  calculate  the  epicadmium 
flux  or  cadmium-covered  detector  response.  Equation  (3-^^)  inay  be  rewritten 
as 


2.(0)     ^  ^ 


e  ^J   C.   +  e  ^J'  S.^ 
jk  jk^ 


%^vV 


(5.1^) 


The  argument  1   has  been  added  to  K  as  a  reminder  that  \  is  defined  as  in 
Eq.  (5.68)  except  that  x  replaces  t  .   Insertion  of  Eqs.  (5-1^)  into  Eq. 

5 

(5.18)  yields 
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2b      2(0)   2 

"0 


3     k 


Mj^(z,t) 


e  ^-^   C,  +  e  ^-^   S., 


;K'\^  •    ^5-15) 


The  response  of  the  l/v  detector  covered  with  cadmium  can  be  taken  to  be 
proportional  to 


u 


epi  —■ ' 


r2E     pcd 


U 


cd 


q(r,T)  e  '   du  , 


(5.16) 


where  u   ,   is   the   lethargy  of  the   cadmiiira  cutoff.      Insertion  of  Eq.    (5-15/ 
into  Eq.    (5 -16)    gives 


,         ,        ,       \     1  (2i    -   l)ny     ^s^Q^      Ti 


i 


x^ 


!„.,    C,    +   I' .,    S., 
^jk     jk         ijk     jk 


;K^\^ 


J     k 


where 


(5.17) 


red  "^-"^        /p 

I^        =      /        p(u^w)    Mj^(z_,u)    e        -^      e^'      du 


(5.18) 


and 
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""cd  -Bl^u 


du   .  (5.19) 


As  before^  experimental  measurements  of  the  quantity  represented  by  Eq. 
(5.17)  taken  as  a  function  of  position  and  frequency  may  be  operated  on 
with  the  operators  given  in  Eqs.  (5. 12)  and  (5.I3).   Let  the  result  of 
operating  with  Eq.  (5-12)  on  Eq.  (5-17)  be  called  F„ .,  ,  which  is 

23(0)     ^      7^,     (2j  -  Drfx 

^^jk=  z^ToT  2  ^  ^ijk^^°^ 2S ^^V^k)  •       (5.20) 


The  nonsymmetric  part  may  be  called  F'   and  is  obtained  by  operating  with 
Eq.  (5.15)  on  Eq.  (5-17),  giving 

2g(0)     ^       7^     jrtx^ 

^m  =  ^;roy  2  i  H^^  j=  ^^^  ^-  ^^v^k^   •  (5.21) 

Several  things  should  be  noted  at  this  point.   First,  all  the  inform- 
ation obtained  thus  far  could  have  been  obtained  from  experimental  data 
taken  at  only  one  z  value.   Secondly,  each  term  obtained  experimentally 
except  for  the  homogeneous  term  has  contained  a  weighted  summation  over  the 
thermal  neutron  flux  at  the  surface  of  the  foils.   Alternatively,  by  com- 
bining the  thermal  constant,  7^,  with  the  flux,  these  summations  may  be 
considered  to  include  the  net  current  into  each  of  the  foils. 

Assvime  now  that  the  experimental  data  and  follo\d.ng  analysis  have 
been  performed  for  M  values  of  z,  different  from  the  locations  of  the  M 
foils,  of  course.   These  M  values  of  z  should  be  chosen  so  as  to  be  fairly 
close  to  the  M  foil  locations  in  order  to  maximize  the  information  about 
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the  foils  in  the  M  sets  of  measurements.   Now  assume  that  this  set  of 
information  is  available  for  the  symmetric  absorption  term,  the  second 
term  in  Eq_.  (5.5).  A  set  of  M  inhomogeneous  equations  results,  with  the  M 
unknowns  being  the  net  current  into  each  of  the  M  foils  or,  if  7  is  known, 
being  the  surface  flux  at  each  of  the  M  foils. 

Wow  that  the  net  c-urrent  into  each  foil  is  known,  this  information 
may  be  inserted  into  some  of  the  other  terms  to  obtain  additional  informa- 
tion.  For  instance,  having  the  net  currents  and  experimental  values  of 
F.  .,  at  M  points,  the  integrals  1„  .,    defined  in  Eq.  (5'l8)  may  be  obtained 
by  solving  the  set  of  M  equations  given  by  Eq.  (5. 20).   Alternatively,  if 
the  E   's  are  assumed  to  be  known,  the  integrals  over  the  slowing-down 
kernel  in  Eq.  (5.5)  may  be  obtained  for  each  foil. 


APPENniX  A 

THE  COMPLEX  AEITHMETIC  VERSION  OF  GIW 

GIM  is  a  FORTRAN  subroutine  developed  by  Burrus  et  al.  (12)  for 
the  purpose  of  obtaining  working  inverses  to  singiolar  matrices  or  to  matrices 
which  appear  singular  to  a  digital  computer.   It  was  modified  to  handle 
complex  arithmetic  and  was  used  in  the  programs  employed  for  calculations 
in  this  dissertation  because,  at  one  point,  it  appeared  that  the  matrices 
in  these  problems  were  very  ill  conditioned.  Although  this  turned  out  not 
to  be  the  case,  it  was  left  in  because  it  takes  up  little  more  storage  than 
a  stajidard  matrix  inverter  and  gives  slightly  better  answers  in  practically 
all  sittiations. 

Theoretically,  the  inverse  of  a  matrix  A,  A~^,  has  the  property  that 
a"-"-  A  =  I,  where  I  is  the  identity  matrix.   This  property  never  is  true 

precisely,  when  an  inverse  is  obtained  on  a  digital  computer  with  finite 

W 
accuracy.   GINV  operates  so  as  to  find  the  working  inverse  A  ,  defined  as 

/     W 
that  matrix  for  which  the  root  mean  square  value  of  (I  -  A  Aj  is  a  minimum. 

The  theory  of  operation  will  not  be  given  here  since  it  is  not 
pertinent  to  this  dissertation  and  should  be  available  soon  in  the 
literature  (ij) . 

Pages  90-92  contain  the  F0RTRAU  listing  of  the  program  as  modified 
for  complex  arithmetic.   Sufficient  comment  cards  are  included  to  enable 
a  competent  programmer  to  modify  the  program  for  ordinary  F0RTRA1T  arithmetic. 
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A  few  comments  may  be  appropriate  on  the  use  of  the  program.   First, 
when  TAU  is  set  to  zero,  the  routine  performs  virtually  the  same  operations 
as  other  matrix  inverters  and  should  give  identical  results  on  well- 
conditioned  matrices.   It  may  also  be  noted  that  if  the  problem  to  be  solved 
is  to  obtain  x  from  Ax  =  b,  where  A  is  an  m  x  n  matrix  (with  n  <  m)  and  b  is 
a  given  m  component  vector,  GHW  automatically  obtains  the  least-squares 
solution  to  the  set  of  n  equations.   It  is  obviously  much  more  versatile 
than  most  matrix  routines,  although  its  fixll  versatility  was  not  needed 
in  this  work. 
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SUaROUTINE  olNV(A,MR,MR,rNC.rAU,UfATtKP) 
C 

C  UPON  LNTKY    A  =  A  MATRIX  WITH  NR  RuWS  Ai^>ID  NC  CULUMNS 

C  AFTER  THL  RIjTURN  THE  ORIGINAL  A  IS  UESTKOYEU  ANU  THE 

C  ARRAY  A  CCNIAINS  THE  TRANSPOSE  OF  T,  THE  TRANSFORMATION 

C  HATRIX  (X  =  r»B)  WHICH  SOLVES  THE  PAIR  OF  EQUATIONS 
C 

C  A»X  =  B       AND  TAU  I»X  =  0 
C 

C  IN  THE  LEAST  SWUARE  SENSE. 
C 

C  MR  =  1ST  DIM.  NO.  OF  ARRAY  A  IN  THE  CALLING  PROGRAMS 

C  TAU  IS  A  NON-NEGATIVE  CONSTANT  WHICH  CONTROLS  THE  ERROR 

C  PROPAGAIIOW  OF  THE  TRANSFORMATION  MATRIX.    IF  TAU  =  0. 

C  THEN  THE  RESULTING  TRANSFORMATION  MATRIX  IS  THE 

C  ORDINARY  LEAST  SUUARES  TRANSFORMATION  MATRIX. 

C  (THE  INVERSE  IF  NR  ^  NO 

C  ATEMP  A.^O  U  ARE  USED  FOR  WORKING  SPACE  BY  THE  ALGORITHM 

C  AND  DO  NOr  NECESSARILY  CONTAIN  ANY  RELEVANT  NUMBERS  AT 

C  THE  CONCLUSION.   ATEMP  MUST  BE  DIMENSIONED  AT  LEAST  NC 

C  BY  THE  CALLING  PROGRAMS.   U  MUST  BE  DIMENSIONED  AT 

C  LEAST  NC»NC  BY  THE  CALLING  PROGRAMS. 

C  NO.  OF  MULTIPLICATIONS  =  NC»*2  1^/2  NR  +  Z/i    NC ) 
C 

I  DIMENSION  A(3U,30),  U(30,30),  ATEMP(30),  DUT(l), 

1   Dur2(l),  TAUSQd),  OUM(l) 
C 

C  THE  DIM.  ST.  AS  SENT  BY  ROSS  BURRUS  FOR  FORTRAN  IV  WAS 

C  DIMENSION  A(MR,NC),U(NC,NC) .ATEMPINC) 
C 

C  NOTE  THAT  IN  COMPLEX  ARITHMETIC  VERSION,  IT  HAS  BEEN 

C  ASSUMED  THAT  MR  IS  ALSO  THE  2N0  DIM.  NO.  OF  A  AS  WELL 

C  AS  THE  1ST. 

rAUSQ(l)  =  rAU»»2 

TAUSg(2)  =  0.0 

C  PLACE  UNIT  MATRIX  IN  U 

DO   5   I  -  l,NC 

DO   4   J  =  1,NC 

I  4          U(I,J)  =  (0.0,0.0) 

I  5               U(I,I)  ^(i. 0,0.0) 
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C      URTHOGGNALI/.L  COMBINED  MATRIX  (A  ABUVt  U)  BY  GRAMM- 

C      SCHMiUT-HILBfcRT  METHOD  WITH  FIRST  NK  ROWS  WEIGHTED 

C      WITH  I  AND  THE  OTHER  NC  ROWS  WEIGHTED  WITH  l/TAU.   THEN 

C      REORfHOGCNALIZE  TO  LESSEN  RGUNOUFF  ERROR. 

C 

DO   20   I  =  1,NC 

IMK  =  I  +  MR 
II  =  I  -  I 
IF   (II)   2,11,2 
2      DC)   10   LL  =  l,.i 

DO   10   J  =  1,11 
JMR  =  J+MR 
I  DOT  =10.0,0.0) 

I  D0T2  =(0.0,0.0) 

DO   3   K  =  1,NR 

DUM( 1)  =  A(K,J) 
DUM(2)  =  -A(K,JMR) 
I    3  .         DOT  =  A(K,  I)»DUM  *■    DOT 

C 
C 

DO   6   K  =  1,J 

DUM( 1)  =  U(K,J) 
DUM(2)  -  -U(K,JMR) 
i    6  D0T2  -    U(K,I)»DUM  +  00T2 

I  DOT  =  DOT  ♦  D0T2»TAUSQ 

DO   8   K  =  l,J 
I    8  U(K,I)  =  U(K,I)  -  DOT»U(K,J) 

DO   10   K  =  1,NK 
I        10  A(K,I)  =  A(K,I)  -  DDT*A(K,J) 

C 

C      NORMALIZE  THE  COLUMN   I  OF  THE  COMBINED  MATRIX 
C 

I   11      DOT  =(0.0,0.0) 
I  DO T2  =(0.0,0.0) 

DO   12   K  =  1,NR 
12  DUI  =  DOT  +  A(K,I)»»2  +  A(K,IMR)»»2 

DO   14   K  =  1,1 
14  0012  =  D0T2  +  U(K,I)»»2  *■    U{K,IMK)»»2 

DOr  =  DOT  *  00T2»TAUSQ 
DOT  =  SURTF(DOT) 
00  17   K  =  1 ,  I 
I   17  U(K,I)  =  U(K, I )/DOT 

DO   19   K  =  1,NR 
I   19  A(K,I)  =  A(K,I )/DOT 

20  CONTINUE 
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U      CALCULATION  UF  THE  IRANSPOSE  UF  THE  TKANSFOKMAT I  ON 

C  MATRIX    TCTRANS.)  =  A  »  U(TR.} 

C 

UO  bO       I  =  l.NR 

UU   45   J  =  l.NC 
I  ATEi'^P(J)  =(0.0,0.0) 

UO   45   K  =  J,NC 
KMR  =  K+MR 
OUMd)  ^  A(  I.K) 
DUM(2)  =  -A{ I, KMR) 
I  ATEMPCJ)  =  UUM«U(J,K)  ♦  ATEMP(J) 

45  CONTINUE 

00   50   J  =  l.NC 
i  A(I,J)  =  ATlMP(J) 

50  CONTlNUt 

RETURN 
END 


APPENDIX  B 


DIGITAL  COMPUTER  TECHNIQUES  FOR  THE  CALCULATION  OF  THE 
ERROR  FUNCTION  WITH  A  COMPLEX  ARGUMENT 


The  calculation  of  the  error  function  vdth  a  complex  argument  is  a 
problem  which  arose  some  time  ago  in  heat  flow  theory  (l^) ,  Two  infinite 
series  approximations  (accurate  to  1  part  in  lO-""®)  are  derived  in  a  paper 
by  Salzer  (15.)  and  two  rather  simple  formulas^  valid  for  large  arguments, 
are  given  by  Gautschi  (l6) .  Two  of  the  above  techniques  have  been  incor- 
porated into  computer  subroutines  that  are  believed  to  give  accuracies  of 
at  least  1  part  in  10®  over  virtually  the  entire  region  of  the  complex 
plane  which  can  be  handled  by  a  56 -bit  digital  computer,  such  as  the 
International  Business  Machine  (IBM)  709- 

Since  Salzer  covers  the  derivation  of  the  infinite  series  approxi- 
mations quite  thoroughly,  only  a  brief  introduction  to  the  development, 
to  illustrate  the  complexity  of  the  problem,  will  be  given.   The  ordinary 
error  function  is  usually  defined  as 

x 
erf(x)  =  --  /   e"'^  dt   .  (B.l) 

^   o 

When  the  complex  argument  is  introduced,  the  integration  must  proceed  along 
some  line  in  the  complex  plane.   If  the  integration  goes  from  the  point 
(0,0)  to  the  point  (X,Y),  one  such  possible  path  is  from  0  to  X  along  the 
real  axis,  then  vertically  from  (X,0)  to  (X,Y) .   This  results  in 


95 


9^ 
X+iY  X  Y 

(B.2) 
The  first  integral  may  be  recognized  as  erf(x),  and  further  manipulation 

with  the  second  integral  gives 


2e-^'  r  rV  ...  o.,. .,_.  r.y^ 


;rf  (X  +  iY)  =  erf  (x)  + ^   /  e"'   sin  2Xy  dy  +  i  /  e"'   cos  2Xy  dy 

(B.5) 


The  two  integrals  in  Eq.  (B.5)  must  be  performed  numerically.   The  very 
clever  technique  described  by  Salzer  results  in  infinite  series  approxi- 
mations to  the  two  functions  in  the  integrals,  each  term  of  which  can  be 
integrated  to  give  infinite  series  approximations  to  the  two  integrals, 
which  are  accurate  to  1  part  in  lO"*"^.   The  result  is 


-X2 
erf(X  +  iY)  =  erf(x)  +  |-^  [(l  -  cos  2XY)  +  i  sin  2XY] 


.X2  °°    -rf/k 

) [F  (X,Y)  +  i  g^(X,Y)]  +  e(X,Y)  ,  {B.h) 

""        ^,  n2  +  UX^   ^  " 

n=l 


where 


F  (X,Y)  =  2X  -  2X  cosh  nY  cos  2XY  +  n  sinh  nY  sin  2XY  (B.5) 


n 


and 


g  (X,Y)  =  2X  cosh  nY  sin  2XY  +  n  sinh  nY  cos  2XY  ,  (B.6) 
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and  the  error  term  is  approximately 


e(X,y) I  z   10'^®  |erf(X  +  iY) 


(B.7) 


The  two  formulas  given  by  Gautschi  are 


•w(z)  =  iz 


O.U6l313^       0.09999216      0.00288389^4- 


>-  z 


0.1901635   z^  -  1.78^+^927   z^  -  5-5253^57  ^ 


where 


+  e(z)   , 


(B.8) 


z  =  X  +  ly 


(B.9) 


wi 


(z)  =   e 


1.2i     ,t-,. 


■x/Jt 


=  e    erfc(-iz)   , 


(B.IO) 


;rfc(z)  =  1  -  erf(z)   , 


(B.ll) 


€(z) I  <  2  X  10"^  for  X  >  3-9  or  y  >  3  , 


(B.12) 


and 


■w(z)  =  iz 


0.312^2^2    _^   0.03176536 


z^  -  0.2752551  z^  -  2.72^7^5 


+  Ti(z)  , 


(B.I3) 


with 


l(z)  I  <  10"®  for  X  >  6  or  y  >  6   . 


(B.lU) 


Equations  (B.^),  (B.8),  (B.IO),  and  (B.I3)  are  all  valid  over  the  entire 
complex  plane^  but  Eq.  (B.IO)  must  be  handled  carefully  to  obtain  erf(z) 
in  terms  of  w(z).   Using  the  property 
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w(z)  =  ^r^  ,  (B.15) 

where  the  bar  indicates  the  complement,  and  Eq.  (B.ll),  Eq.  (B.IO)  may  be 
manip\iLated  to  obtain  an  expression  valid  for  x  >  0  and  y  >  0: 


erf(z 


)  =  1  -  e^   (cos  2xy  -  i  sin  2xy)   w(y  +  ix)   .  (B.16) 


Using  Eq.  (B.l6)  to  perform  calculations  in  the  first  quadrant  of  the 
complex  plane,  the  symmetry  properties  of  the  error  function  can  be  used  to 
reach  the  other  three  quadrants.  Defining 


z  =  x  +  iy  , 


z  =  X  -  ly 


X  +  iy 


(B.IT) 
(B.18) 

(B.19) 


the  error  function  has  the  following  properties: 

erf(-z)  =  -  erf(z)   ,  (B.20) 

erf(z)  =  erf(z)   ,  (B.21) 

erf(z)  =  erf(z)   .  (B.22) 

Figure  15  shows  the  various  regions  of  the  complex  plane  when  dif- 
ferent techniques  were  used  to  calculate  the  error  function.   In  Region  I, 
defined  by  |x|  <9.58  and  |y|  <9.58,  Eqs.  (B.h) ,    (B.5),  axid  (B.6)  were 
used.   In  order  to  obtain  an  accuracy  of  1  part  in  10^,  it  was  found  nec- 
essary to  use  approximately  9  +  2|y|  terms  in  Eq.  (B.i+).   In  Region  II, 
defined  by  y^  -  x^  <  -  88.028,  the  error  function  is  (l,0)  for  x  >  0  and 
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Fig.  15.   Regions  of  Complex  Plane  Requiring  Different  Techniques  for 
Calculating  the  Error  Fimction. 
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is  (-1^0)  for  X  <  0,  to  the  numerical  limits  of  an  IBM  709  o^   similar 
computer.   In  Region  III,  defined  by  y^  -  x^  >  88.027,  the  imaginary  part 
(and  sometimes  also  the  real  part)  of  the  error  function  exceeds  the  capacity 
of  an  IBM  709  and  therefore  cannot  be  calculated.   Region  IV  is  defined  by 
the  remaining  four  portions  of  the  complex  plane,  and  in  this  area  Eqs. 
(B.I5)  and  (B.I6)  were  used. 

The  following  two  pages  give  FORTRAN  listings  of  the  two  computer 
subroutines.   SUBR0UTrNE  CMERPN  is  entered  initially.   If  |x|  <  9.58  and 
|y|  <  9-58,  it  uses  Eqs.  (B.U),  (b.5),  and  (B.6)  to  perform  the  calculation 
and  then  returns.   This  calculation  takes  a  considerable  amount  of  machine 
time,  the  measured  running  times  for  the  IBM  709  being  approximately  O.I5 
sec  for  11  terras  and  0.5^  sec  for  28  terms.   If  x  or  y  is  outside  Region  I 
of  Fig.  B.l,  CMERFW  transfers  control  to  SUBROUTINE  ERFGR6.   ERFGR6  gives 
an  error  return  if  it  finds  that  the  point  is  in  Region  III,  and  uses  Eqs. 
(B.I5)  and  (B.I6)  and  the  symmetry  properties  for  Regions  II  and  IV. 

It  may  be  noted  that  ERFGR6  is  a  rather  straightforward  conversion 
of  the  appropriate  equations  to  F0RTRM  statements,  but  it  was  found  nec- 
essary to  reformulate  Eqs.  (B.U),  (B.5),  and  (B.6)  slightly  for  CMERFN 
in  order  to  cover  the  f\ill  range  of  y  values.   The  form  used  was  obtained 
by  expressing  the  hyperbolic  functions  in  Eqs.  (B.5)  and  (B.6)  in  terms  of 
exponentials,  then  combining  these  exponentials  with  the  factor  exp(-n^/U). 
Another  noteworthy  feature  of  this  subroutine  is  that  it  uses  four  terms  of 
the  Maclaurin  series  for  the  terms  (l  -  cos  2xy) /x  and  (sin  2xy) /x  when 
|xy|  <  0.1,  rather  than  using  the  FORTRAN  SIN  and  Cl^S  functions.   This 
feature  was  found  to  be  necessary  to  ensirre  8-digit  accuracy  for  the 
complete  region. 
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SUaRUUTlNE  CMERFN  ( X, Y , REERFN.CMER , N ) 
C   FOR  ABSF(X)  AND  ABSF(Y)  LESS  THAN  9.38,  THIS  RQUTI-NE  USES 
C     UP  TU  28  TERMS  OF  A  ShRKS  TO  CALCULATE  THE  ERF(X*IY). 
C   OUTSIDE  THIS  RANGE,  IT  CALLS  ERFGR6. 

N=0 

IF  (ABSF(X)-9.38)  10,10,15 
10  IF  (A8SF{Y)-9.38)  25,25,15 
15  CALL  ERFGR6(X,Y,REERFN,CMER) 

RETURN 
25  CONS  =  .636619772 
C     CONS  =  2/PI 

CONS  =  CONS/EXPF{X»X) 

ERFOFX  =  LRFF(X) 

TWOX  =  2.0»X 

C0S2XY  =  CUSF  ( TWOX»Y) 

SXN2XY  =  SINF  (TWOX»Y) 

IF  (ABSF{X*Y)-0.1)  50,50,75 
50  XY2  =  X»X»Y»Y 

A  =  X»Y»Y»{.5-XY2/6.»( l.-XY2/7.5»(l.-XY2/14.)) ) 

8  =  Y»( .5-XY2/3.»(l.-XY2/5.»(l.-XY2/l0.5)) ) 

GO  ru  100 

75  A  =  (1.0-CUS2XY)/2./TWOX 

B  =  SIN2XY/2,/TW0X 
100  N  =  N+1 

IF  (N-30)  125,125,175 
1/5  REERFN  =    ERFOFX  +  CUNS»A 

CMER  =  CONS»B 

WRITE  OUTPUT  TAPE  6, 150 , N, REERFN ,CMtR, AERR , BERR, X,Y 
150  FORMAT  (I10,6E18.9) 

RETURN 
125  FN  =  FLOATF(N) 

FNY  =  FN»Y 

DUMY=1./{FN»*2+TW0X»«2) 
255  EXPP  =  tXPF(FN*(Y-FN/4. } )/2. 

EXPN  =  tXPF(-FN»(Y+FN/4.) )/2. 

ADUM=£XPP»(FN»SIN2XY-TWGX«C0S2XY) 

1  -EXPN»(FN»SIN2XY+TW0X«CUS2XY) 

2  +lW0X»£XPF(-FN»FN/4.) 
QDUM=EXPP»(FN»C0S2XY+TW0X»SIN2XY) 

1       -EXPN»(FN*C0S2XY-TW0X»SIN2XY) 
275  ADUM  =  DUMY*ADUM 

BDUM  =  UUMY»BDUM 

AERR  =  ADUM/A 

BERR  =  ROUM/a 

A  =  A4-AUUM 

B  =  B+BDUM 
300  IF  (ABSF(AtRR)  ♦■  ABSF(BtRR)  -  l.E-8)  200,200,100 
200  CONTINUE 

REERFN=  ERFOFX  ♦  CONS»A 

CMER   =  CUNS*B 

RETURN 

END 
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c 
c 
c 
I 


SU3  PvOUT  I  \L-  f.  RFGK6  (  X  ,  Y  ,  Z  I  ,  Z2  ) 
IF  Y»«2-X*»2  IS  GRLATER  THA;^  e<;.027,  AiM 
blVtN.   If-  Y»»2-X»»2  IS  LLSS  THAN  -88 
RETURN  SIbNF(1.0,X)  AS  Li    ^Hl>    ZERO  AS 
DIKlNSION  Z{  1)  ,o(  1  ) 
IF{Y»Y-X«X-88.027)  b  ,  i) ,  3 
3  WRiTL  UUTPUT  TAHE  6,4, X,Y 
A  FtJKMAl  (2'^H   CMEKFN  UNABLE  10  HANDLE  X 
2     9\\  UR  Y  =fc20.9) 

REIU^^ 
5  A  =  (0. '3124242,0.) 
B  =  {0.2^b2  5l3l,(;.  ) 
(0.05176536,0. 
(2.  724745,0.  ) 
SlGNFl  1.0,X) 
:  0  ,  V  ) 


ERROR  RETURM  IS 
028,  IT  SHOULD 
IZ. 


^E20.9, 


) 


C  = 

D  = 

SX  = 

SY  =  SIUNF(  1 

X  =  AI5SF{X) 

Y  =  ABSF(Y) 

/{  1)  =  Y 

Z(2)  =  X 

G  =  1*7 

Z  =  A/{b-D)  +  C/(G-U) 

G(  1  )  =  -X 

G  (  2  )  =  Y 

Z  =    G*Z 

Z12)  =  -Z(2) 

LX  =  lX?F( Y*Y-X»X) 

^  .:X«CGSF(  2.»X»Y) 
--    -LX»SI;jF(  2.»X«Y) 
,0,0.0) 
G«Z 


G(  1  ) 
G(  2  ) 
A  =  (  +  1, 

Z  =  ^  - 


Z  I  ^  Z  (  1  )  *  S  X 
IZ    =  Z{2)*^SY 

KETUKJ 
END 


APPENDIX  C 


GALMDI'S  THERMAL  CONSTAJM  FOR  A 
SIWUSOIDALLY  VARYING  FLUX 


The  thermal  constant^  y ,    defined  by  Feinberg  and  Galanin  as  the 
ratio  of  the  net  current  into  a  small  slxig  or  foil  to  the  neutron  flux  on 
its  surface,  is  a  real  constant  in  the  steady-state  case.   It  is  fairly  easy 
to  show  that  7  becomes  complex  when  the  flux  is  varying  sinusoidally.   If 
diffusion  theory  is  applicable  in  the  slug,  it  is  evident  that  the  factor 
1/L^  in  the  usual  diffusion  equation  must  be  replaced  by  1/l^  +  i  oj/vD, 

Using  diffusion  theory,  Galanin  (17)  obtains  the  following  formulas. 
For  an  infinite  cylinder  of  radius  p, 

where  D  and  L  are  the  diffusion  coefficient  and  diffusion  length  in  the 
foil,  and  Iq  and  Ii  are  modified  Bessel  functions  of  the  first  kind.   For 
an  infinite  slab  of  thickness  2d, 


L  coth  f - 


It  should  be  noted  that  the  above  equation  includes  a  factor  of  2  which  was 
left  out  of  Galanin 's  form-ula  (56.2). 
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The  above  formulas  are  applicable  for  the  present  work  simply  by 
making  the  substitution  noted  above  for  L,    if  one  is  willing  to  make  the 
assumptions  of  diffusion  theory  and  infinite  dimensions  of  the  fuel  slugs 
or  foils.  More  elaborate  calculations  are  certainly  feasible  althoiigh 
somewhat  involved. 


APPENDIX  D 

PHYSICAL  CONSTANTS  USED  IN  SAMPLE  CALCULATIONS 

In  order  to  permit  reproduction  of  the  sample  calculations  of 
Chapter  IV,  the  physical  constants  which  were  used  are  given  below: 


Parameter 

Graphite 

Heavy  Water 
[39 -^i   EfeO) 

P(g 

0.95 

0.999 

T 
s 

I.U7  X  10"* 

5.08  X  10"^ 

T 
S 

36U. 

125. 

z 

a 

U.09  X  lO"'* 

1.57  X  10'* 

D 

0.986 

0.825 

p  (t  )  and  T  are  defined  in  Eqs.  (U.5)  and  (U.5).  t  is  the  Fermi  age 

So  " 

to  thermal,  in  cm^,  as  defined  by  Eq.  (5.II).   Z  and  D  are  defined  in 
Eq.  (1.1). 

The  thermal  constant,  y,  was  taken  to  be  7.8252  for  the  cadmium 
foils  and  O.5227  for  the  uranium  foils.  r\  for  uranium  was  taken  to  be 
2.07. 
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